Actual source code: da2.c

      1: #define PETSCDM_DLL
      2:
     3:  #include src/dm/da/daimpl.h

      7: /*@C
      8:       DAGetElements - Gets an array containing the indices (in local coordinates)
      9:                  of all the local elements

     11:     Not Collective

     13:    Input Parameter:
     14: .     da - the DA object

     16:    Output Parameters:
     17: +     n - number of local elements
     18: -     e - the indices of the elements vertices

     20:    Level: intermediate

     22: .seealso: DAElementType, DASetElementType(), DARestoreElements()
     23: @*/
     24: PetscErrorCode PETSCDM_DLLEXPORT DAGetElements(DA da,PetscInt *n,const PetscInt *e[])
     25: {
     29:   (da->ops->getelements)(da,n,e);
     30:   return(0);
     31: }

     35: /*@C
     36:       DARestoreElements - Returns an array containing the indices (in local coordinates)
     37:                  of all the local elements obtained with DAGetElements()

     39:     Not Collective

     41:    Input Parameter:
     42: +     da - the DA object
     43: .     n - number of local elements
     44: -     e - the indices of the elements vertices

     46:    Level: intermediate

     48: .seealso: DAElementType, DASetElementType(), DAGetElements()
     49: @*/
     50: PetscErrorCode PETSCDM_DLLEXPORT DARestoreElements(DA da,PetscInt *n,const PetscInt *e[])
     51: {
     55:   if (da->ops->restoreelements) {
     56:     (da->ops->restoreelements)(da,n,e);
     57:   }
     58:   return(0);
     59: }

     63: PetscErrorCode PETSCDM_DLLEXPORT DAGetOwnershipRange(DA da,PetscInt **lx,PetscInt **ly,PetscInt **lz)
     64: {
     67:   if (lx) *lx = da->lx;
     68:   if (ly) *ly = da->ly;
     69:   if (lz) *lz = da->lz;
     70:   return(0);
     71: }

     75: PetscErrorCode DAView_2d(DA da,PetscViewer viewer)
     76: {
     78:   PetscMPIInt    rank;
     79:   PetscTruth     iascii,isdraw;

     82:   MPI_Comm_rank(da->comm,&rank);

     84:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
     85:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);
     86:   if (iascii) {
     87:     PetscViewerASCIISynchronizedPrintf(viewer,"Processor [%d] M %D N %D m %D n %D w %D s %D\n",rank,da->M,
     88:                              da->N,da->m,da->n,da->w,da->s);
     89:     PetscViewerASCIISynchronizedPrintf(viewer,"X range of indices: %D %D, Y range of indices: %D %D\n",da->xs,da->xe,da->ys,da->ye);
     90:     PetscViewerFlush(viewer);
     91:   } else if (isdraw) {
     92:     PetscDraw       draw;
     93:     double     ymin = -1*da->s-1,ymax = da->N+da->s;
     94:     double     xmin = -1*da->s-1,xmax = da->M+da->s;
     95:     double     x,y;
     96:     PetscInt   base,*idx;
     97:     char       node[10];
     98:     PetscTruth isnull;
     99:
    100:     PetscViewerDrawGetDraw(viewer,0,&draw);
    101:     PetscDrawIsNull(draw,&isnull); if (isnull) return(0);
    102:     if (!da->coordinates) {
    103:       PetscDrawSetCoordinates(draw,xmin,ymin,xmax,ymax);
    104:     }
    105:     PetscDrawSynchronizedClear(draw);

    107:     /* first processor draw all node lines */
    108:     if (!rank) {
    109:       ymin = 0.0; ymax = da->N - 1;
    110:       for (xmin=0; xmin<da->M; xmin++) {
    111:         PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_BLACK);
    112:       }
    113:       xmin = 0.0; xmax = da->M - 1;
    114:       for (ymin=0; ymin<da->N; ymin++) {
    115:         PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_BLACK);
    116:       }
    117:     }
    118:     PetscDrawSynchronizedFlush(draw);
    119:     PetscDrawPause(draw);

    121:     /* draw my box */
    122:     ymin = da->ys; ymax = da->ye - 1; xmin = da->xs/da->w;
    123:     xmax =(da->xe-1)/da->w;
    124:     PetscDrawLine(draw,xmin,ymin,xmax,ymin,PETSC_DRAW_RED);
    125:     PetscDrawLine(draw,xmin,ymin,xmin,ymax,PETSC_DRAW_RED);
    126:     PetscDrawLine(draw,xmin,ymax,xmax,ymax,PETSC_DRAW_RED);
    127:     PetscDrawLine(draw,xmax,ymin,xmax,ymax,PETSC_DRAW_RED);

    129:     /* put in numbers */
    130:     base = (da->base)/da->w;
    131:     for (y=ymin; y<=ymax; y++) {
    132:       for (x=xmin; x<=xmax; x++) {
    133:         sprintf(node,"%d",(int)base++);
    134:         PetscDrawString(draw,x,y,PETSC_DRAW_BLACK,node);
    135:       }
    136:     }

    138:     PetscDrawSynchronizedFlush(draw);
    139:     PetscDrawPause(draw);
    140:     /* overlay ghost numbers, useful for error checking */
    141:     /* put in numbers */

    143:     base = 0; idx = da->idx;
    144:     ymin = da->Ys; ymax = da->Ye; xmin = da->Xs; xmax = da->Xe;
    145:     for (y=ymin; y<ymax; y++) {
    146:       for (x=xmin; x<xmax; x++) {
    147:         if ((base % da->w) == 0) {
    148:           sprintf(node,"%d",(int)(idx[base]/da->w));
    149:           PetscDrawString(draw,x/da->w,y,PETSC_DRAW_BLUE,node);
    150:         }
    151:         base++;
    152:       }
    153:     }
    154:     PetscDrawSynchronizedFlush(draw);
    155:     PetscDrawPause(draw);
    156:   } else {
    157:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for DA2d",((PetscObject)viewer)->type_name);
    158:   }
    159:   return(0);
    160: }

    164: PetscErrorCode DAPublish_Petsc(PetscObject obj)
    165: {
    167:   return(0);
    168: }

    172: PetscErrorCode DAGetElements_2d_P1(DA da,PetscInt *n,const PetscInt *e[])
    173: {
    175:   PetscInt       i,j,cnt,xs,xe = da->xe,ys,ye = da->ye,Xs = da->Xs, Xe = da->Xe, Ys = da->Ys;

    178:   if (!da->e) {
    179:     if (da->xs == Xs) xs = da->xs; else xs = da->xs - 1;
    180:     if (da->ys == Ys) ys = da->ys; else ys = da->ys - 1;
    181:     da->ne = 2*(xe - xs - 1)*(ye - ys - 1);
    182:     PetscMalloc((1 + 3*da->ne)*sizeof(PetscInt),&da->e);
    183:     cnt    = 0;
    184:     for (j=ys; j<ye-1; j++) {
    185:       for (i=xs; i<xe-1; i++) {
    186:         da->e[cnt]   = i - Xs + (j - Ys)*(Xe - Xs);
    187:         da->e[cnt+1] = i - Xs + 1 + (j - Ys)*(Xe - Xs);
    188:         da->e[cnt+2] = i - Xs + (j - Ys + 1)*(Xe - Xs);

    190:         da->e[cnt+3] = i - Xs + 1 + (j - Ys + 1)*(Xe - Xs);
    191:         da->e[cnt+4] = i - Xs + (j - Ys + 1)*(Xe - Xs);
    192:         da->e[cnt+5] = i - Xs + 1 + (j - Ys)*(Xe - Xs);
    193:         cnt += 6;
    194:       }
    195:     }
    196:   }
    197:   *n = da->ne;
    198:   *e = da->e;
    199:   return(0);
    200: }

    205: /*@C
    206:    DACreate2d -  Creates an object that will manage the communication of  two-dimensional
    207:    regular array data that is distributed across some processors.

    209:    Collective on MPI_Comm

    211:    Input Parameters:
    212: +  comm - MPI communicator
    213: .  wrap - type of periodicity should the array have.
    214:          Use one of DA_NONPERIODIC, DA_XPERIODIC, DA_YPERIODIC, or DA_XYPERIODIC.
    215: .  stencil_type - stencil type.  Use either DA_STENCIL_BOX or DA_STENCIL_STAR.
    216: .  M,N - global dimension in each direction of the array (use -M and or -N to indicate that it may be set to a different value
    217:             from the command line with -da_grid_x <M> -da_grid_y <N>)
    218: .  m,n - corresponding number of processors in each dimension
    219:          (or PETSC_DECIDE to have calculated)
    220: .  dof - number of degrees of freedom per node
    221: .  s - stencil width
    222: -  lx, ly - arrays containing the number of nodes in each cell along
    223:            the x and y coordinates, or PETSC_NULL. If non-null, these
    224:            must be of length as m and n, and the corresponding
    225:            m and n cannot be PETSC_DECIDE. The sum of the lx[] entries
    226:            must be M, and the sum of the ly[] entries must be N.

    228:    Output Parameter:
    229: .  inra - the resulting distributed array object

    231:    Options Database Key:
    232: +  -da_view - Calls DAView() at the conclusion of DACreate2d()
    233: .  -da_grid_x <nx> - number of grid points in x direction, if M < 0
    234: .  -da_grid_y <ny> - number of grid points in y direction, if N < 0
    235: .  -da_processors_x <nx> - number of processors in x direction
    236: .  -da_processors_y <ny> - number of processors in y direction
    237: .  -da_refine_x - refinement ratio in x direction
    238: -  -da_refine_y - refinement ratio in y direction

    240:    Level: beginner

    242:    Notes:
    243:    The stencil type DA_STENCIL_STAR with width 1 corresponds to the
    244:    standard 5-pt stencil, while DA_STENCIL_BOX with width 1 denotes
    245:    the standard 9-pt stencil.

    247:    The array data itself is NOT stored in the DA, it is stored in Vec objects;
    248:    The appropriate vector objects can be obtained with calls to DACreateGlobalVector()
    249:    and DACreateLocalVector() and calls to VecDuplicate() if more are needed.

    251: .keywords: distributed array, create, two-dimensional

    253: .seealso: DADestroy(), DAView(), DACreate1d(), DACreate3d(), DAGlobalToLocalBegin(), DAGetRefinementFactor(),
    254:           DAGlobalToLocalEnd(), DALocalToGlobal(), DALocalToLocalBegin(), DALocalToLocalEnd(), DASetRefinementFactor(),
    255:           DAGetInfo(), DACreateGlobalVector(), DACreateLocalVector(), DACreateNaturalVector(), DALoad(), DAView()

    257: @*/
    258: PetscErrorCode PETSCDM_DLLEXPORT DACreate2d(MPI_Comm comm,DAPeriodicType wrap,DAStencilType stencil_type,
    259:                           PetscInt M,PetscInt N,PetscInt m,PetscInt n,PetscInt dof,PetscInt s,PetscInt *lx,PetscInt *ly,DA *inra)
    260: {
    262:   PetscMPIInt    rank,size;
    263:   PetscInt       xs,xe,ys,ye,x,y,Xs,Xe,Ys,Ye,start,end;
    264:   PetscInt       up,down,left,i,n0,n1,n2,n3,n5,n6,n7,n8,*idx,nn;
    265:   PetscInt       xbase,*bases,*ldims,j,x_t,y_t,s_t,base,count;
    266:   PetscInt       s_x,s_y; /* s proportionalized to w */
    267:   PetscInt       *flx = 0,*fly = 0;
    268:   PetscInt       sn0 = 0,sn2 = 0,sn6 = 0,sn8 = 0,refine_x = 2, refine_y = 2,tM = M,tN = N;
    269:   PetscTruth     flg1;
    270:   DA             da;
    271:   Vec            local,global;
    272:   VecScatter     ltog,gtol;
    273:   IS             to,from;

    277:   *inra = 0;
    278: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
    279:   DMInitializePackage(PETSC_NULL);
    280: #endif

    282:   if (dof < 1) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Must have 1 or more degrees of freedom per node: %D",dof);
    283:   if (s < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Stencil width cannot be negative: %D",s);

    285:   PetscOptionsBegin(comm,PETSC_NULL,"2d DA Options","DA");
    286:     if (M < 0){
    287:       tM = -M;
    288:       PetscOptionsInt("-da_grid_x","Number of grid points in x direction","DACreate2d",tM,&tM,PETSC_NULL);
    289:     }
    290:     if (N < 0){
    291:       tN = -N;
    292:       PetscOptionsInt("-da_grid_y","Number of grid points in y direction","DACreate2d",tN,&tN,PETSC_NULL);
    293:     }
    294:     PetscOptionsInt("-da_processors_x","Number of processors in x direction","DACreate2d",m,&m,PETSC_NULL);
    295:     PetscOptionsInt("-da_processors_y","Number of processors in y direction","DACreate2d",n,&n,PETSC_NULL);
    296:     PetscOptionsInt("-da_refine_x","Refinement ratio in x direction","DASetRefinementFactor",refine_x,&refine_x,PETSC_NULL);
    297:     PetscOptionsInt("-da_refine_y","Refinement ratio in y direction","DASetRefinementFactor",refine_y,&refine_y,PETSC_NULL);
    298:   PetscOptionsEnd();
    299:   M = tM; N = tN;

    301:   PetscHeaderCreate(da,_p_DA,struct _DAOps,DA_COOKIE,0,"DA",comm,DADestroy,DAView);
    302:   da->bops->publish           = DAPublish_Petsc;
    303:   da->ops->createglobalvector = DACreateGlobalVector;
    304:   da->ops->getinterpolation   = DAGetInterpolation;
    305:   da->ops->getcoloring        = DAGetColoring;
    306:   da->ops->getmatrix          = DAGetMatrix;
    307:   da->ops->refine             = DARefine;
    308:   da->ops->getinjection       = DAGetInjection;
    309:   da->ops->getelements        = DAGetElements_2d_P1;
    310:   da->elementtype             = DA_ELEMENT_P1;

    312:   PetscLogObjectMemory(da,sizeof(struct _p_DA));
    313:   da->dim        = 2;
    314:   da->interptype = DA_Q1;
    315:   da->refine_x   = refine_x;
    316:   da->refine_y   = refine_y;
    317:   PetscMalloc(dof*sizeof(char*),&da->fieldname);
    318:   PetscMemzero(da->fieldname,dof*sizeof(char*));

    320:   MPI_Comm_size(comm,&size);
    321:   MPI_Comm_rank(comm,&rank);

    323:   if (m != PETSC_DECIDE) {
    324:     if (m < 1) {SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in X direction: %D",m);}
    325:     else if (m > size) {SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in X direction: %D %d",m,size);}
    326:   }
    327:   if (n != PETSC_DECIDE) {
    328:     if (n < 1) {SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Non-positive number of processors in Y direction: %D",n);}
    329:     else if (n > size) {SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Too many processors in Y direction: %D %d",n,size);}
    330:   }

    332:   if (m == PETSC_DECIDE || n == PETSC_DECIDE) {
    333:     /* try for squarish distribution */
    334:     /* This should use MPI_Dims_create instead */
    335:     m = (PetscInt)(0.5 + sqrt(((double)M)*((double)size)/((double)N)));
    336:     if (!m) m = 1;
    337:     while (m > 0) {
    338:       n = size/m;
    339:       if (m*n == size) break;
    340:       m--;
    341:     }
    342:     if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;}
    343:     if (m*n != size) SETERRQ(PETSC_ERR_PLIB,"Internally Created Bad Partition");
    344:   } else if (m*n != size) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Given Bad partition");

    346:   if (M < m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Partition in x direction is too fine! %D %D",M,m);
    347:   if (N < n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Partition in y direction is too fine! %D %D",N,n);

    349:   /*
    350:      We should create an MPI Cartesian topology here, with reorder
    351:      set to true.  That would create a NEW communicator that we would
    352:      need to use for operations on this distributed array
    353:   */

    355:   /*
    356:      Determine locally owned region
    357:      xs is the first local node number, x is the number of local nodes
    358:   */
    359:   if (!lx) { /* user sets distribution */
    360:     PetscMalloc(m*sizeof(PetscInt),&lx);
    361:     flx = lx;
    362:     for (i=0; i<m; i++) {
    363:       lx[i] = M/m + ((M % m) > i);
    364:     }
    365:   }
    366:   x  = lx[rank % m];
    367:   xs = 0;
    368:   for (i=0; i<(rank % m); i++) {
    369:     xs += lx[i];
    370:   }
    371: #if defined(PETSC_USE_DEBUGGING)
    372:   left = xs;
    373:   for (i=(rank % m); i<m; i++) {
    374:     left += lx[i];
    375:   }
    376:   if (left != M) {
    377:     SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Sum of lx across processors not equal to M: %D %D",left,M);
    378:   }
    379: #endif

    381:   /*
    382:      Determine locally owned region
    383:      ys is the first local node number, y is the number of local nodes
    384:   */
    385:   if (!ly) { /* user sets distribution */
    386:     PetscMalloc(n*sizeof(PetscInt),&ly);
    387:     fly  = ly;
    388:     for (i=0; i<n; i++) {
    389:       ly[i] = N/n + ((N % n) > i);
    390:     }
    391:   }
    392:   y  = ly[rank/m];
    393:   ys = 0;
    394:   for (i=0; i<(rank/m); i++) {
    395:     ys += ly[i];
    396:   }
    397: #if defined(PETSC_USE_DEBUGGING)
    398:   left = ys;
    399:   for (i=(rank/m); i<n; i++) {
    400:     left += ly[i];
    401:   }
    402:   if (left != N) {
    403:     SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Sum of ly across processors not equal to N: %D %D",left,N);
    404:   }
    405: #endif

    407:   xe = xs + x;
    408:   ye = ys + y;

    410:   /* determine ghost region */
    411:   /* Assume No Periodicity */
    412:   if (xs-s > 0) Xs = xs - s; else Xs = 0;
    413:   if (ys-s > 0) Ys = ys - s; else Ys = 0;
    414:   if (xe+s <= M) Xe = xe + s; else Xe = M;
    415:   if (ye+s <= N) Ye = ye + s; else Ye = N;

    417:   /* X Periodic */
    418:   if (DAXPeriodic(wrap)){
    419:     Xs = xs - s;
    420:     Xe = xe + s;
    421:   }

    423:   /* Y Periodic */
    424:   if (DAYPeriodic(wrap)){
    425:     Ys = ys - s;
    426:     Ye = ye + s;
    427:   }

    429:   /* Resize all X parameters to reflect w */
    430:   x   *= dof;
    431:   xs  *= dof;
    432:   xe  *= dof;
    433:   Xs  *= dof;
    434:   Xe  *= dof;
    435:   s_x = s*dof;
    436:   s_y = s;

    438:   /* determine starting point of each processor */
    439:   nn    = x*y;
    440:   PetscMalloc((2*size+1)*sizeof(PetscInt),&bases);
    441:   ldims = bases+size+1;
    442:   MPI_Allgather(&nn,1,MPIU_INT,ldims,1,MPIU_INT,comm);
    443:   bases[0] = 0;
    444:   for (i=1; i<=size; i++) {
    445:     bases[i] = ldims[i-1];
    446:   }
    447:   for (i=1; i<=size; i++) {
    448:     bases[i] += bases[i-1];
    449:   }

    451:   /* allocate the base parallel and sequential vectors */
    452:   da->Nlocal = x*y;
    453:   VecCreateMPIWithArray(comm,da->Nlocal,PETSC_DECIDE,0,&global);
    454:   VecSetBlockSize(global,dof);
    455:   da->nlocal = (Xe-Xs)*(Ye-Ys);
    456:   VecCreateSeqWithArray(PETSC_COMM_SELF,da->nlocal,0,&local);
    457:   VecSetBlockSize(local,dof);

    460:   /* generate appropriate vector scatters */
    461:   /* local to global inserts non-ghost point region into global */
    462:   VecGetOwnershipRange(global,&start,&end);
    463:   ISCreateStride(comm,x*y,start,1,&to);

    465:   left  = xs - Xs; down  = ys - Ys; up    = down + y;
    466:   PetscMalloc(x*(up - down)*sizeof(PetscInt),&idx);
    467:   count = 0;
    468:   for (i=down; i<up; i++) {
    469:     for (j=0; j<x/dof; j++) {
    470:       idx[count++] = left + i*(Xe-Xs) + j*dof;
    471:     }
    472:   }
    473:   ISCreateBlock(comm,dof,count,idx,&from);
    474:   PetscFree(idx);

    476:   VecScatterCreate(local,from,global,to,&ltog);
    477:   PetscLogObjectParent(da,to);
    478:   PetscLogObjectParent(da,from);
    479:   PetscLogObjectParent(da,ltog);
    480:   ISDestroy(from);
    481:   ISDestroy(to);

    483:   /* global to local must include ghost points */
    484:   if (stencil_type == DA_STENCIL_BOX) {
    485:     ISCreateStride(comm,(Xe-Xs)*(Ye-Ys),0,1,&to);
    486:   } else {
    487:     /* must drop into cross shape region */
    488:     /*       ---------|
    489:             |  top    |
    490:          |---         ---|
    491:          |   middle      |
    492:          |               |
    493:          ----         ----
    494:             | bottom  |
    495:             -----------
    496:         Xs xs        xe  Xe */
    497:     /* bottom */
    498:     left  = xs - Xs; down = ys - Ys; up    = down + y;
    499:     count = down*(xe-xs) + (up-down)*(Xe-Xs) + (Ye-Ys-up)*(xe-xs);
    500:     PetscMalloc(count*sizeof(PetscInt)/dof,&idx);
    501:     count = 0;
    502:     for (i=0; i<down; i++) {
    503:       for (j=0; j<xe-xs; j += dof) {
    504:         idx[count++] = left + i*(Xe-Xs) + j;
    505:       }
    506:     }
    507:     /* middle */
    508:     for (i=down; i<up; i++) {
    509:       for (j=0; j<Xe-Xs; j += dof) {
    510:         idx[count++] = i*(Xe-Xs) + j;
    511:       }
    512:     }
    513:     /* top */
    514:     for (i=up; i<Ye-Ys; i++) {
    515:       for (j=0; j<xe-xs; j += dof) {
    516:         idx[count++] = left + i*(Xe-Xs) + j;
    517:       }
    518:     }
    519:     ISCreateBlock(comm,dof,count,idx,&to);
    520:     PetscFree(idx);
    521:   }

    524:   /* determine who lies on each side of us stored in    n6 n7 n8
    525:                                                         n3    n5
    526:                                                         n0 n1 n2
    527:   */

    529:   /* Assume the Non-Periodic Case */
    530:   n1 = rank - m;
    531:   if (rank % m) {
    532:     n0 = n1 - 1;
    533:   } else {
    534:     n0 = -1;
    535:   }
    536:   if ((rank+1) % m) {
    537:     n2 = n1 + 1;
    538:     n5 = rank + 1;
    539:     n8 = rank + m + 1; if (n8 >= m*n) n8 = -1;
    540:   } else {
    541:     n2 = -1; n5 = -1; n8 = -1;
    542:   }
    543:   if (rank % m) {
    544:     n3 = rank - 1;
    545:     n6 = n3 + m; if (n6 >= m*n) n6 = -1;
    546:   } else {
    547:     n3 = -1; n6 = -1;
    548:   }
    549:   n7 = rank + m; if (n7 >= m*n) n7 = -1;

    552:   /* Modify for Periodic Cases */
    553:   if (wrap == DA_YPERIODIC) {  /* Handle Top and Bottom Sides */
    554:     if (n1 < 0) n1 = rank + m * (n-1);
    555:     if (n7 < 0) n7 = rank - m * (n-1);
    556:     if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
    557:     if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1;
    558:     if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
    559:     if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1;
    560:   } else if (wrap == DA_XPERIODIC) { /* Handle Left and Right Sides */
    561:     if (n3 < 0) n3 = rank + (m-1);
    562:     if (n5 < 0) n5 = rank - (m-1);
    563:     if ((n1 >= 0) && (n0 < 0)) n0 = rank-1;
    564:     if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1;
    565:     if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1;
    566:     if ((n7 >= 0) && (n8 < 0)) n8 = rank+1;
    567:   } else if (wrap == DA_XYPERIODIC) {

    569:     /* Handle all four corners */
    570:     if ((n6 < 0) && (n7 < 0) && (n3 < 0)) n6 = m-1;
    571:     if ((n8 < 0) && (n7 < 0) && (n5 < 0)) n8 = 0;
    572:     if ((n2 < 0) && (n5 < 0) && (n1 < 0)) n2 = size-m;
    573:     if ((n0 < 0) && (n3 < 0) && (n1 < 0)) n0 = size-1;

    575:     /* Handle Top and Bottom Sides */
    576:     if (n1 < 0) n1 = rank + m * (n-1);
    577:     if (n7 < 0) n7 = rank - m * (n-1);
    578:     if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1;
    579:     if ((n3 >= 0) && (n6 < 0)) n6 = (rank%m)-1;
    580:     if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1;
    581:     if ((n5 >= 0) && (n8 < 0)) n8 = (rank%m)+1;

    583:     /* Handle Left and Right Sides */
    584:     if (n3 < 0) n3 = rank + (m-1);
    585:     if (n5 < 0) n5 = rank - (m-1);
    586:     if ((n1 >= 0) && (n0 < 0)) n0 = rank-1;
    587:     if ((n1 >= 0) && (n2 < 0)) n2 = rank-2*m+1;
    588:     if ((n7 >= 0) && (n6 < 0)) n6 = rank+2*m-1;
    589:     if ((n7 >= 0) && (n8 < 0)) n8 = rank+1;
    590:   }

    592:   if (stencil_type == DA_STENCIL_STAR) {
    593:     /* save corner processor numbers */
    594:     sn0 = n0; sn2 = n2; sn6 = n6; sn8 = n8;
    595:     n0 = n2 = n6 = n8 = -1;
    596:   }

    598:   PetscMalloc((x+2*s_x)*(y+2*s_y)*sizeof(PetscInt),&idx);
    599:   PetscLogObjectMemory(da,(x+2*s_x)*(y+2*s_y)*sizeof(PetscInt));
    600:   nn = 0;

    602:   xbase = bases[rank];
    603:   for (i=1; i<=s_y; i++) {
    604:     if (n0 >= 0) { /* left below */
    605:       x_t = lx[n0 % m]*dof;
    606:       y_t = ly[(n0/m)];
    607:       s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x;
    608:       for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    609:     }
    610:     if (n1 >= 0) { /* directly below */
    611:       x_t = x;
    612:       y_t = ly[(n1/m)];
    613:       s_t = bases[n1] + x_t*y_t - (s_y+1-i)*x_t;
    614:       for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    615:     }
    616:     if (n2 >= 0) { /* right below */
    617:       x_t = lx[n2 % m]*dof;
    618:       y_t = ly[(n2/m)];
    619:       s_t = bases[n2] + x_t*y_t - (s_y+1-i)*x_t;
    620:       for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    621:     }
    622:   }

    624:   for (i=0; i<y; i++) {
    625:     if (n3 >= 0) { /* directly left */
    626:       x_t = lx[n3 % m]*dof;
    627:       /* y_t = y; */
    628:       s_t = bases[n3] + (i+1)*x_t - s_x;
    629:       for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    630:     }

    632:     for (j=0; j<x; j++) { idx[nn++] = xbase++; } /* interior */

    634:     if (n5 >= 0) { /* directly right */
    635:       x_t = lx[n5 % m]*dof;
    636:       /* y_t = y; */
    637:       s_t = bases[n5] + (i)*x_t;
    638:       for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    639:     }
    640:   }

    642:   for (i=1; i<=s_y; i++) {
    643:     if (n6 >= 0) { /* left above */
    644:       x_t = lx[n6 % m]*dof;
    645:       /* y_t = ly[(n6/m)]; */
    646:       s_t = bases[n6] + (i)*x_t - s_x;
    647:       for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    648:     }
    649:     if (n7 >= 0) { /* directly above */
    650:       x_t = x;
    651:       /* y_t = ly[(n7/m)]; */
    652:       s_t = bases[n7] + (i-1)*x_t;
    653:       for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    654:     }
    655:     if (n8 >= 0) { /* right above */
    656:       x_t = lx[n8 % m]*dof;
    657:       /* y_t = ly[(n8/m)]; */
    658:       s_t = bases[n8] + (i-1)*x_t;
    659:       for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    660:     }
    661:   }

    663:   base = bases[rank];
    664:   {
    665:     PetscInt nnn = nn/dof,*iidx;
    666:     PetscMalloc(nnn*sizeof(PetscInt),&iidx);
    667:     for (i=0; i<nnn; i++) {
    668:       iidx[i] = idx[dof*i];
    669:     }
    670:     ISCreateBlock(comm,dof,nnn,iidx,&from);
    671:     PetscFree(iidx);
    672:   }
    673:   VecScatterCreate(global,from,local,to,&gtol);
    674:   PetscLogObjectParent(da,to);
    675:   PetscLogObjectParent(da,from);
    676:   PetscLogObjectParent(da,gtol);
    677:   ISDestroy(to);
    678:   ISDestroy(from);

    680:   if (stencil_type == DA_STENCIL_STAR) {
    681:     /*
    682:         Recompute the local to global mappings, this time keeping the
    683:       information about the cross corner processor numbers.
    684:     */
    685:     n0 = sn0; n2 = sn2; n6 = sn6; n8 = sn8;
    686:     nn = 0;
    687:     xbase = bases[rank];
    688:     for (i=1; i<=s_y; i++) {
    689:       if (n0 >= 0) { /* left below */
    690:         x_t = lx[n0 % m]*dof;
    691:         y_t = ly[(n0/m)];
    692:         s_t = bases[n0] + x_t*y_t - (s_y-i)*x_t - s_x;
    693:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    694:       }
    695:       if (n1 >= 0) { /* directly below */
    696:         x_t = x;
    697:         y_t = ly[(n1/m)];
    698:         s_t = bases[n1] + x_t*y_t - (s_y+1-i)*x_t;
    699:         for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    700:       }
    701:       if (n2 >= 0) { /* right below */
    702:         x_t = lx[n2 % m]*dof;
    703:         y_t = ly[(n2/m)];
    704:         s_t = bases[n2] + x_t*y_t - (s_y+1-i)*x_t;
    705:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    706:       }
    707:     }

    709:     for (i=0; i<y; i++) {
    710:       if (n3 >= 0) { /* directly left */
    711:         x_t = lx[n3 % m]*dof;
    712:         /* y_t = y; */
    713:         s_t = bases[n3] + (i+1)*x_t - s_x;
    714:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    715:       }

    717:       for (j=0; j<x; j++) { idx[nn++] = xbase++; } /* interior */

    719:       if (n5 >= 0) { /* directly right */
    720:         x_t = lx[n5 % m]*dof;
    721:         /* y_t = y; */
    722:         s_t = bases[n5] + (i)*x_t;
    723:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    724:       }
    725:     }

    727:     for (i=1; i<=s_y; i++) {
    728:       if (n6 >= 0) { /* left above */
    729:         x_t = lx[n6 % m]*dof;
    730:         /* y_t = ly[(n6/m)]; */
    731:         s_t = bases[n6] + (i)*x_t - s_x;
    732:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    733:       }
    734:       if (n7 >= 0) { /* directly above */
    735:         x_t = x;
    736:         /* y_t = ly[(n7/m)]; */
    737:         s_t = bases[n7] + (i-1)*x_t;
    738:         for (j=0; j<x_t; j++) { idx[nn++] = s_t++;}
    739:       }
    740:       if (n8 >= 0) { /* right above */
    741:         x_t = lx[n8 % m]*dof;
    742:         /* y_t = ly[(n8/m)]; */
    743:         s_t = bases[n8] + (i-1)*x_t;
    744:         for (j=0; j<s_x; j++) { idx[nn++] = s_t++;}
    745:       }
    746:     }
    747:   }
    748:   PetscFree(bases);

    750:   da->M  = M;  da->N  = N;  da->m  = m;  da->n  = n;  da->w = dof;  da->s = s;
    751:   da->xs = xs; da->xe = xe; da->ys = ys; da->ye = ye; da->zs = 0; da->ze = 1;
    752:   da->Xs = Xs; da->Xe = Xe; da->Ys = Ys; da->Ye = Ye; da->Zs = 0; da->Ze = 1;
    753:   da->P  = 1;  da->p  = 1;

    755:   VecDestroy(local);
    756:   VecDestroy(global);

    758:   da->gtol         = gtol;
    759:   da->ltog         = ltog;
    760:   da->idx          = idx;
    761:   da->Nl           = nn;
    762:   da->base         = base;
    763:   da->wrap         = wrap;
    764:   da->ops->view    = DAView_2d;
    765:   da->stencil_type = stencil_type;

    767:   /*
    768:      Set the local to global ordering in the global vector, this allows use
    769:      of VecSetValuesLocal().
    770:   */
    771:   ISLocalToGlobalMappingCreateNC(comm,nn,idx,&da->ltogmap);
    772:   ISLocalToGlobalMappingBlock(da->ltogmap,da->w,&da->ltogmapb);
    773:   PetscLogObjectParent(da,da->ltogmap);

    775:   *inra = da;

    777:   da->ltol = PETSC_NULL;
    778:   da->ao   = PETSC_NULL;

    781:   if (!flx) {
    782:     PetscMalloc(m*sizeof(PetscInt),&flx);
    783:     PetscMemcpy(flx,lx,m*sizeof(PetscInt));
    784:   }
    785:   if (!fly) {
    786:     PetscMalloc(n*sizeof(PetscInt),&fly);
    787:     PetscMemcpy(fly,ly,n*sizeof(PetscInt));
    788:   }
    789:   da->lx = flx;
    790:   da->ly = fly;

    792:   PetscOptionsHasName(PETSC_NULL,"-da_view",&flg1);
    793:   if (flg1) {DAView(da,PETSC_VIEWER_STDOUT_(da->comm));}
    794:   PetscOptionsHasName(PETSC_NULL,"-da_view_draw",&flg1);
    795:   if (flg1) {DAView(da,PETSC_VIEWER_DRAW_(da->comm));}
    796:   PetscOptionsHasName(PETSC_NULL,"-help",&flg1);
    797:   if (flg1) {DAPrintHelp(da);}

    799:   PetscPublishAll(da);
    800:   return(0);
    801: }

    805: /*@
    806:    DAPrintHelp - Prints command line options for DA.

    808:    Collective on DA

    810:    Input Parameters:
    811: .  da - the distributed array

    813:    Level: intermediate

    815: .seealso: DACreate1d(), DACreate2d(), DACreate3d()

    817: .keywords: DA, help

    819: @*/
    820: PetscErrorCode PETSCDM_DLLEXPORT DAPrintHelp(DA da)
    821: {
    822:   static PetscTruth called = PETSC_FALSE;
    823:   MPI_Comm          comm;
    824:   PetscErrorCode    ierr;

    829:   comm = da->comm;
    830:   if (!called) {
    831:     (*PetscHelpPrintf)(comm,"General Distributed Array (DA) options:\n");
    832:     (*PetscHelpPrintf)(comm,"  -da_view: print DA distribution to screen\n");
    833:     (*PetscHelpPrintf)(comm,"  -da_view_draw: display DA in window\n");
    834:     called = PETSC_TRUE;
    835:   }
    836:   return(0);
    837: }

    841: /*@C
    842:    DARefine - Creates a new distributed array that is a refinement of a given
    843:    distributed array.

    845:    Collective on DA

    847:    Input Parameter:
    848: +  da - initial distributed array
    849: -  comm - communicator to contain refined DA, must be either same as the da communicator or include the
    850:           da communicator and be 2, 4, or 8 times larger. Currently ignored

    852:    Output Parameter:
    853: .  daref - refined distributed array

    855:    Level: advanced

    857:    Note:
    858:    Currently, refinement consists of just doubling the number of grid spaces
    859:    in each dimension of the DA.

    861: .keywords:  distributed array, refine

    863: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy()
    864: @*/
    865: PetscErrorCode PETSCDM_DLLEXPORT DARefine(DA da,MPI_Comm comm,DA *daref)
    866: {
    868:   PetscInt       M,N,P;
    869:   DA             da2;

    875:   if (DAXPeriodic(da->wrap) || da->interptype == DA_Q0){
    876:     M = da->refine_x*da->M;
    877:   } else {
    878:     M = 1 + da->refine_x*(da->M - 1);
    879:   }
    880:   if (DAYPeriodic(da->wrap) || da->interptype == DA_Q0){
    881:     N = da->refine_y*da->N;
    882:   } else {
    883:     N = 1 + da->refine_y*(da->N - 1);
    884:   }
    885:   if (DAZPeriodic(da->wrap) || da->interptype == DA_Q0){
    886:     P = da->refine_z*da->P;
    887:   } else {
    888:     P = 1 + da->refine_z*(da->P - 1);
    889:   }
    890:   if (da->dim == 1) {
    891:     DACreate1d(da->comm,da->wrap,M,da->w,da->s,PETSC_NULL,&da2);
    892:   } else if (da->dim == 2) {
    893:     DACreate2d(da->comm,da->wrap,da->stencil_type,M,N,da->m,da->n,da->w,da->s,PETSC_NULL,PETSC_NULL,&da2);
    894:   } else if (da->dim == 3) {
    895:     DACreate3d(da->comm,da->wrap,da->stencil_type,M,N,P,da->m,da->n,da->p,da->w,da->s,0,0,0,&da2);
    896:   }
    897:   /* allow overloaded (user replaced) operations to be inherited by refinement clones */
    898:   da2->ops->getmatrix        = da->ops->getmatrix;
    899:   da2->ops->getinterpolation = da->ops->getinterpolation;
    900:   da2->ops->getcoloring      = da->ops->getcoloring;
    901:   da2->interptype            = da->interptype;
    902:
    903:   /* copy fill information if given */
    904:   if (da->dfill) {
    905:     PetscMalloc((da->dfill[da->w]+da->w+1)*sizeof(PetscInt),&da2->dfill);
    906:     PetscMemcpy(da2->dfill,da->dfill,(da->dfill[da->w]+da->w+1)*sizeof(PetscInt));
    907:   }
    908:   if (da->ofill) {
    909:     PetscMalloc((da->ofill[da->w]+da->w+1)*sizeof(PetscInt),&da2->ofill);
    910:     PetscMemcpy(da2->ofill,da->ofill,(da->ofill[da->w]+da->w+1)*sizeof(PetscInt));
    911:   }
    912:   /* copy the refine information */
    913:   da2->refine_x = da->refine_x;
    914:   da2->refine_y = da->refine_y;
    915:   da2->refine_z = da->refine_z;
    916:   *daref = da2;
    917:   return(0);
    918: }

    920: /*@C
    921:      DASetRefinementFactor - Set the ratios that the DA grid is refined

    923:     Collective on DA

    925:   Input Parameters:
    926: +    da - the DA object
    927: .    refine_x - ratio of fine grid to coarse in x direction (2 by default)
    928: .    refine_y - ratio of fine grid to coarse in y direction (2 by default)
    929: -    refine_z - ratio of fine grid to coarse in z direction (2 by default)

    931:   Options Database:
    932: +  -da_refine_x - refinement ratio in x direction
    933: .  -da_refine_y - refinement ratio in y direction
    934: -  -da_refine_y - refinement ratio in z direction

    936:   Level: intermediate

    938:     Notes: Pass PETSC_IGNORE to leave a value unchanged

    940: .seealso: DARefine(), DAGetRefinementFactor()
    941: @*/
    942: PetscErrorCode PETSCDM_DLLEXPORT DASetRefinementFactor(DA da, PetscInt refine_x, PetscInt refine_y,PetscInt refine_z)
    943: {
    945:   if (refine_x > 0) da->refine_x = refine_x;
    946:   if (refine_y > 0) da->refine_y = refine_y;
    947:   if (refine_z > 0) da->refine_z = refine_z;
    948:   return(0);
    949: }

    951: /*@C
    952:      DAGetRefinementFactor - Gets the ratios that the DA grid is refined

    954:     Not Collective

    956:   Input Parameter:
    957: .    da - the DA object

    959:   Output Parameters:
    960: +    refine_x - ratio of fine grid to coarse in x direction (2 by default)
    961: .    refine_y - ratio of fine grid to coarse in y direction (2 by default)
    962: -    refine_z - ratio of fine grid to coarse in z direction (2 by default)

    964:   Level: intermediate

    966:     Notes: Pass PETSC_NULL for values you do not need

    968: .seealso: DARefine(), DASetRefinementFactor()
    969: @*/
    970: PetscErrorCode PETSCDM_DLLEXPORT DAGetRefinementFactor(DA da, PetscInt *refine_x, PetscInt *refine_y,PetscInt *refine_z)
    971: {
    973:   if (refine_x) *refine_x = da->refine_x;
    974:   if (refine_y) *refine_y = da->refine_y;
    975:   if (refine_z) *refine_z = da->refine_z;
    976:   return(0);
    977: }

    979: /*@C
    980:      DASetGetMatrix - Sets the routine used by the DA to allocate a matrix.

    982:     Collective on DA

    984:   Input Parameters:
    985: +    da - the DA object
    986: -    f - the function that allocates the matrix for that specific DA

    988:   Level: developer

    990:    Notes: See DASetBlockFills() that provides a simple way to provide the nonzero structure for
    991:        the diagonal and off-diagonal blocks of the matrix

    993: .seealso: DAGetMatrix(), DASetBlockFills()
    994: @*/
    995: PetscErrorCode PETSCDM_DLLEXPORT DASetGetMatrix(DA da,PetscErrorCode (*f)(DA, MatType,Mat*))
    996: {
    998:   da->ops->getmatrix = f;
    999:   return(0);
    1000: }

    1002: /*
    1003:       M is number of grid points
    1004:       m is number of processors

    1006: */
    1009: PetscErrorCode PETSCDM_DLLEXPORT DASplitComm2d(MPI_Comm comm,PetscInt M,PetscInt N,PetscInt sw,MPI_Comm *outcomm)
    1010: {
    1012:   PetscInt       m,n = 0,x = 0,y = 0;
    1013:   PetscMPIInt    size,csize,rank;

    1016:   MPI_Comm_size(comm,&size);
    1017:   MPI_Comm_rank(comm,&rank);

    1019:   csize = 4*size;
    1020:   do {
    1021:     if (csize % 4) SETERRQ4(PETSC_ERR_ARG_INCOMP,"Cannot split communicator of size %d tried %d %D %D",size,csize,x,y);
    1022:     csize   = csize/4;
    1023:
    1024:     m = (PetscInt)(0.5 + sqrt(((double)M)*((double)csize)/((double)N)));
    1025:     if (!m) m = 1;
    1026:     while (m > 0) {
    1027:       n = csize/m;
    1028:       if (m*n == csize) break;
    1029:       m--;
    1030:     }
    1031:     if (M > N && m < n) {PetscInt _m = m; m = n; n = _m;}

    1033:     x = M/m + ((M % m) > ((csize-1) % m));
    1034:     y = (N + (csize-1)/m)/n;
    1035:   } while ((x < 4 || y < 4) && csize > 1);
    1036:   if (size != csize) {
    1037:     MPI_Group    entire_group,sub_group;
    1038:     PetscMPIInt  i,*groupies;

    1040:     MPI_Comm_group(comm,&entire_group);
    1041:     PetscMalloc(csize*sizeof(PetscInt),&groupies);
    1042:     for (i=0; i<csize; i++) {
    1043:       groupies[i] = (rank/csize)*csize + i;
    1044:     }
    1045:     MPI_Group_incl(entire_group,csize,groupies,&sub_group);
    1046:     PetscFree(groupies);
    1047:     MPI_Comm_create(comm,sub_group,outcomm);
    1048:     MPI_Group_free(&entire_group);
    1049:     MPI_Group_free(&sub_group);
    1050:     PetscLogInfo((0,"DASplitComm2d:Creating redundant coarse problems of size %d\n",csize));
    1051:   } else {
    1052:     *outcomm = comm;
    1053:   }
    1054:   return(0);
    1055: }

    1059: /*@C
    1060:        DASetLocalFunction - Caches in a DA a local function.

    1062:    Collective on DA

    1064:    Input Parameter:
    1065: +  da - initial distributed array
    1066: -  lf - the local function

    1068:    Level: intermediate

    1070:    Notes: The routine SNESDAFormFunction() uses this the cached function to evaluate the user provided function.

    1072: .keywords:  distributed array, refine

    1074: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DAGetLocalFunction(), DASetLocalFunctioni()
    1075: @*/
    1076: PetscErrorCode PETSCDM_DLLEXPORT DASetLocalFunction(DA da,DALocalFunction1 lf)
    1077: {
    1080:   da->lf    = lf;
    1081:   return(0);
    1082: }

    1086: /*@C
    1087:        DASetLocalFunctioni - Caches in a DA a local function that evaluates a single component

    1089:    Collective on DA

    1091:    Input Parameter:
    1092: +  da - initial distributed array
    1093: -  lfi - the local function

    1095:    Level: intermediate

    1097: .keywords:  distributed array, refine

    1099: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DAGetLocalFunction(), DASetLocalFunction()
    1100: @*/
    1101: PetscErrorCode PETSCDM_DLLEXPORT DASetLocalFunctioni(DA da,PetscErrorCode (*lfi)(DALocalInfo*,MatStencil*,void*,PetscScalar*,void*))
    1102: {
    1105:   da->lfi = lfi;
    1106:   return(0);
    1107: }

    1112: PetscErrorCode DASetLocalAdicFunction_Private(DA da,DALocalFunction1 ad_lf)
    1113: {
    1116:   da->adic_lf = ad_lf;
    1117:   return(0);
    1118: }

    1120: /*MC
    1121:        DASetLocalAdicFunctioni - Caches in a DA a local functioni computed by ADIC/ADIFOR

    1123:    Collective on DA

    1125:    Synopsis:
    1126:    PetscErrorCode DASetLocalAdicFunctioni(DA da,PetscInt (ad_lf*)(DALocalInfo*,MatStencil*,void*,void*,void*)
    1127:
    1128:    Input Parameter:
    1129: +  da - initial distributed array
    1130: -  ad_lfi - the local function as computed by ADIC/ADIFOR

    1132:    Level: intermediate

    1134: .keywords:  distributed array, refine

    1136: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DAGetLocalFunction(), DASetLocalFunction(),
    1137:           DASetLocalJacobian(), DASetLocalFunctioni()
    1138: M*/

    1142: PetscErrorCode DASetLocalAdicFunctioni_Private(DA da,PetscErrorCode (*ad_lfi)(DALocalInfo*,MatStencil*,void*,void*,void*))
    1143: {
    1146:   da->adic_lfi = ad_lfi;
    1147:   return(0);
    1148: }

    1150: /*MC
    1151:        DASetLocalAdicMFFunctioni - Caches in a DA a local functioni computed by ADIC/ADIFOR

    1153:    Collective on DA

    1155:    Synopsis:
    1156:    PetscErrorCode  DASetLocalAdicFunctioni(DA da,int (ad_lf*)(DALocalInfo*,MatStencil*,void*,void*,void*)
    1157:
    1158:    Input Parameter:
    1159: +  da - initial distributed array
    1160: -  admf_lfi - the local matrix-free function as computed by ADIC/ADIFOR

    1162:    Level: intermediate

    1164: .keywords:  distributed array, refine

    1166: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DAGetLocalFunction(), DASetLocalFunction(),
    1167:           DASetLocalJacobian(), DASetLocalFunctioni()
    1168: M*/

    1172: PetscErrorCode DASetLocalAdicMFFunctioni_Private(DA da,PetscErrorCode (*admf_lfi)(DALocalInfo*,MatStencil*,void*,void*,void*))
    1173: {
    1176:   da->adicmf_lfi = admf_lfi;
    1177:   return(0);
    1178: }

    1180: /*MC
    1181:        DASetLocalAdicMFFunction - Caches in a DA a local function computed by ADIC/ADIFOR

    1183:    Collective on DA

    1185:    Synopsis:
    1186:    PetscErrorCode DASetLocalAdicMFFunction(DA da,DALocalFunction1 ad_lf)
    1187:
    1188:    Input Parameter:
    1189: +  da - initial distributed array
    1190: -  ad_lf - the local function as computed by ADIC/ADIFOR

    1192:    Level: intermediate

    1194: .keywords:  distributed array, refine

    1196: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DAGetLocalFunction(), DASetLocalFunction(),
    1197:           DASetLocalJacobian()
    1198: M*/

    1202: PetscErrorCode DASetLocalAdicMFFunction_Private(DA da,DALocalFunction1 ad_lf)
    1203: {
    1206:   da->adicmf_lf = ad_lf;
    1207:   return(0);
    1208: }

    1210: /*@C
    1211:        DASetLocalJacobian - Caches in a DA a local Jacobian

    1213:    Collective on DA

    1215:
    1216:    Input Parameter:
    1217: +  da - initial distributed array
    1218: -  lj - the local Jacobian

    1220:    Level: intermediate

    1222:    Notes: The routine SNESDAFormFunction() uses this the cached function to evaluate the user provided function.

    1224: .keywords:  distributed array, refine

    1226: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DAGetLocalFunction(), DASetLocalFunction()
    1227: @*/
    1230: PetscErrorCode PETSCDM_DLLEXPORT DASetLocalJacobian(DA da,DALocalFunction1 lj)
    1231: {
    1234:   da->lj    = lj;
    1235:   return(0);
    1236: }

    1240: /*@C
    1241:        DAGetLocalFunction - Gets from a DA a local function and its ADIC/ADIFOR Jacobian

    1243:    Collective on DA

    1245:    Input Parameter:
    1246: .  da - initial distributed array

    1248:    Output Parameters:
    1249: .  lf - the local function

    1251:    Level: intermediate

    1253: .keywords:  distributed array, refine

    1255: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DASetLocalFunction()
    1256: @*/
    1257: PetscErrorCode PETSCDM_DLLEXPORT DAGetLocalFunction(DA da,DALocalFunction1 *lf)
    1258: {
    1261:   if (lf)       *lf = da->lf;
    1262:   return(0);
    1263: }

    1267: /*@
    1268:     DAFormFunction1 - Evaluates a user provided function on each processor that
    1269:         share a DA

    1271:    Input Parameters:
    1272: +    da - the DA that defines the grid
    1273: .    vu - input vector
    1274: .    vfu - output vector
    1275: -    w - any user data

    1277:     Notes: Does NOT do ghost updates on vu upon entry

    1279:     Level: advanced

    1281: .seealso: DAComputeJacobian1WithAdic()

    1283: @*/
    1284: PetscErrorCode PETSCDM_DLLEXPORT DAFormFunction1(DA da,Vec vu,Vec vfu,void *w)
    1285: {
    1287:   void           *u,*fu;
    1288:   DALocalInfo    info;
    1289:

    1292:   DAGetLocalInfo(da,&info);
    1293:   DAVecGetArray(da,vu,&u);
    1294:   DAVecGetArray(da,vfu,&fu);

    1296:   (*da->lf)(&info,u,fu,w);
    1297:   if (PetscExceptionValue(ierr)) {
    1298:     PetscErrorCode pDAVecRestoreArray(da,vu,&u);CHKERRQ(pierr);
    1299:     pDAVecRestoreArray(da,vfu,&fu);CHKERRQ(pierr);
    1300:   }
    1301:

    1303:   DAVecRestoreArray(da,vu,&u);
    1304:   DAVecRestoreArray(da,vfu,&fu);
    1305:   return(0);
    1306: }

    1310: PetscErrorCode PETSCDM_DLLEXPORT DAFormFunctioniTest1(DA da,void *w)
    1311: {
    1312:   Vec            vu,fu,fui;
    1314:   PetscInt       i,n;
    1315:   PetscScalar    *ui,mone = -1.0;
    1316:   PetscRandom    rnd;
    1317:   PetscReal      norm;

    1320:   DAGetLocalVector(da,&vu);
    1321:   PetscRandomCreate(PETSC_COMM_SELF,RANDOM_DEFAULT,&rnd);
    1322:   VecSetRandom(vu,rnd);
    1323:   PetscRandomDestroy(rnd);

    1325:   DAGetGlobalVector(da,&fu);
    1326:   DAGetGlobalVector(da,&fui);
    1327:
    1328:   DAFormFunction1(da,vu,fu,w);

    1330:   VecGetArray(fui,&ui);
    1331:   VecGetLocalSize(fui,&n);
    1332:   for (i=0; i<n; i++) {
    1333:     DAFormFunctioni1(da,i,vu,ui+i,w);
    1334:   }
    1335:   VecRestoreArray(fui,&ui);

    1337:   VecAXPY(fui,mone,fu);
    1338:   VecNorm(fui,NORM_2,&norm);
    1339:   PetscPrintf(da->comm,"Norm of difference in vectors %g\n",norm);
    1340:   VecView(fu,0);
    1341:   VecView(fui,0);

    1343:   DARestoreLocalVector(da,&vu);
    1344:   DARestoreGlobalVector(da,&fu);
    1345:   DARestoreGlobalVector(da,&fui);
    1346:   return(0);
    1347: }

    1351: /*@
    1352:     DAFormFunctioni1 - Evaluates a user provided function

    1354:    Input Parameters:
    1355: +    da - the DA that defines the grid
    1356: .    i - the component of the function we wish to compute (must be local)
    1357: .    vu - input vector
    1358: .    vfu - output value
    1359: -    w - any user data

    1361:     Notes: Does NOT do ghost updates on vu upon entry

    1363:     Level: advanced

    1365: .seealso: DAComputeJacobian1WithAdic()

    1367: @*/
    1368: PetscErrorCode PETSCDM_DLLEXPORT DAFormFunctioni1(DA da,PetscInt i,Vec vu,PetscScalar *vfu,void *w)
    1369: {
    1371:   void           *u;
    1372:   DALocalInfo    info;
    1373:   MatStencil     stencil;
    1374:

    1377:   DAGetLocalInfo(da,&info);
    1378:   DAVecGetArray(da,vu,&u);

    1380:   /* figure out stencil value from i */
    1381:   stencil.c = i % info.dof;
    1382:   stencil.i = (i % (info.xm*info.dof))/info.dof;
    1383:   stencil.j = (i % (info.xm*info.ym*info.dof))/(info.xm*info.dof);
    1384:   stencil.k = i/(info.xm*info.ym*info.dof);

    1386:   (*da->lfi)(&info,&stencil,u,vfu,w);

    1388:   DAVecRestoreArray(da,vu,&u);
    1389:   return(0);
    1390: }

    1392: #if defined(new)
    1395: /*
    1396:   DAGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix where local
    1397:     function lives on a DA

    1399:         y ~= (F(u + ha) - F(u))/h,
    1400:   where F = nonlinear function, as set by SNESSetFunction()
    1401:         u = current iterate
    1402:         h = difference interval
    1403: */
    1404: PetscErrorCode DAGetDiagonal_MFFD(DA da,Vec U,Vec a)
    1405: {
    1406:   PetscScalar    h,*aa,*ww,v;
    1407:   PetscReal      epsilon = PETSC_SQRT_MACHINE_EPSILON,umin = 100.0*PETSC_SQRT_MACHINE_EPSILON;
    1409:   PetscInt       gI,nI;
    1410:   MatStencil     stencil;
    1411:   DALocalInfo    info;
    1412:
    1414:   (*ctx->func)(0,U,a,ctx->funcctx);
    1415:   (*ctx->funcisetbase)(U,ctx->funcctx);

    1417:   VecGetArray(U,&ww);
    1418:   VecGetArray(a,&aa);
    1419:
    1420:   nI = 0;
    1421:     h  = ww[gI];
    1422:     if (h == 0.0) h = 1.0;
    1423: #if !defined(PETSC_USE_COMPLEX)
    1424:     if (h < umin && h >= 0.0)      h = umin;
    1425:     else if (h < 0.0 && h > -umin) h = -umin;
    1426: #else
    1427:     if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0)     h = umin;
    1428:     else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin;
    1429: #endif
    1430:     h     *= epsilon;
    1431:
    1432:     ww[gI += h;
    1433:     (*ctx->funci)(i,w,&v,ctx->funcctx);
    1434:     aa[nI]  = (v - aa[nI])/h;
    1435:     ww[gI] -= h;
    1436:     nI++;
    1437:   }
    1438:   VecRestoreArray(U,&ww);
    1439:   VecRestoreArray(a,&aa);
    1440:   return(0);
    1441: }
    1442: #endif

    1444: #if defined(PETSC_HAVE_ADIC)
    1446: #include "adic/ad_utils.h"

    1451: /*@C
    1452:     DAComputeJacobian1WithAdic - Evaluates a adiC provided Jacobian function on each processor that
    1453:         share a DA

    1455:    Input Parameters:
    1456: +    da - the DA that defines the grid
    1457: .    vu - input vector (ghosted)
    1458: .    J - output matrix
    1459: -    w - any user data

    1461:    Level: advanced

    1463:     Notes: Does NOT do ghost updates on vu upon entry

    1465: .seealso: DAFormFunction1()

    1467: @*/
    1468: PetscErrorCode PETSCDM_DLLEXPORT DAComputeJacobian1WithAdic(DA da,Vec vu,Mat J,void *w)
    1469: {
    1471:   PetscInt       gtdof,tdof;
    1472:   PetscScalar    *ustart;
    1473:   DALocalInfo    info;
    1474:   void           *ad_u,*ad_f,*ad_ustart,*ad_fstart;
    1475:   ISColoring     iscoloring;

    1478:   DAGetLocalInfo(da,&info);

    1480:   PetscADResetIndep();

    1482:   /* get space for derivative objects.  */
    1483:   DAGetAdicArray(da,PETSC_TRUE,(void **)&ad_u,&ad_ustart,&gtdof);
    1484:   DAGetAdicArray(da,PETSC_FALSE,(void **)&ad_f,&ad_fstart,&tdof);
    1485:   VecGetArray(vu,&ustart);
    1486:   DAGetColoring(da,IS_COLORING_GHOSTED,&iscoloring);

    1488:   PetscADSetValueAndColor(ad_ustart,gtdof,iscoloring->colors,ustart);

    1490:   VecRestoreArray(vu,&ustart);
    1491:   ISColoringDestroy(iscoloring);
    1492:   PetscADIncrementTotalGradSize(iscoloring->n);
    1493:   PetscADSetIndepDone();

    1495:   PetscLogEventBegin(DA_LocalADFunction,0,0,0,0);
    1496:   (*da->adic_lf)(&info,ad_u,ad_f,w);
    1497:   PetscLogEventEnd(DA_LocalADFunction,0,0,0,0);

    1499:   /* stick the values into the matrix */
    1500:   MatSetValuesAdic(J,(PetscScalar**)ad_fstart);

    1502:   /* return space for derivative objects.  */
    1503:   DARestoreAdicArray(da,PETSC_TRUE,(void **)&ad_u,&ad_ustart,&gtdof);
    1504:   DARestoreAdicArray(da,PETSC_FALSE,(void **)&ad_f,&ad_fstart,&tdof);
    1505:   return(0);
    1506: }

    1510: /*@C
    1511:     DAMultiplyByJacobian1WithAdic - Applies an ADIC-provided Jacobian function to a vector on
    1512:     each processor that shares a DA.

    1514:     Input Parameters:
    1515: +   da - the DA that defines the grid
    1516: .   vu - Jacobian is computed at this point (ghosted)
    1517: .   v - product is done on this vector (ghosted)
    1518: .   fu - output vector = J(vu)*v (not ghosted)
    1519: -   w - any user data

    1521:     Notes:
    1522:     This routine does NOT do ghost updates on vu upon entry.

    1524:    Level: advanced

    1526: .seealso: DAFormFunction1()

    1528: @*/
    1529: PetscErrorCode PETSCDM_DLLEXPORT DAMultiplyByJacobian1WithAdic(DA da,Vec vu,Vec v,Vec f,void *w)
    1530: {
    1532:   PetscInt       i,gtdof,tdof;
    1533:   PetscScalar    *avu,*av,*af,*ad_vustart,*ad_fstart;
    1534:   DALocalInfo    info;
    1535:   void           *ad_vu,*ad_f;

    1538:   DAGetLocalInfo(da,&info);

    1540:   /* get space for derivative objects.  */
    1541:   DAGetAdicMFArray(da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,&gtdof);
    1542:   DAGetAdicMFArray(da,PETSC_FALSE,(void **)&ad_f,(void**)&ad_fstart,&tdof);

    1544:   /* copy input vector into derivative object */
    1545:   VecGetArray(vu,&avu);
    1546:   VecGetArray(v,&av);
    1547:   for (i=0; i<gtdof; i++) {
    1548:     ad_vustart[2*i]   = avu[i];
    1549:     ad_vustart[2*i+1] = av[i];
    1550:   }
    1551:   VecRestoreArray(vu,&avu);
    1552:   VecRestoreArray(v,&av);

    1554:   PetscADResetIndep();
    1555:   PetscADIncrementTotalGradSize(1);
    1556:   PetscADSetIndepDone();

    1558:   (*da->adicmf_lf)(&info,ad_vu,ad_f,w);

    1560:   /* stick the values into the vector */
    1561:   VecGetArray(f,&af);
    1562:   for (i=0; i<tdof; i++) {
    1563:     af[i] = ad_fstart[2*i+1];
    1564:   }
    1565:   VecRestoreArray(f,&af);

    1567:   /* return space for derivative objects.  */
    1568:   DARestoreAdicMFArray(da,PETSC_TRUE,(void **)&ad_vu,(void**)&ad_vustart,&gtdof);
    1569:   DARestoreAdicMFArray(da,PETSC_FALSE,(void **)&ad_f,(void**)&ad_fstart,&tdof);
    1570:   return(0);
    1571: }
    1572: #endif

    1576: /*@
    1577:     DAComputeJacobian1 - Evaluates a local Jacobian function on each processor that
    1578:         share a DA

    1580:    Input Parameters:
    1581: +    da - the DA that defines the grid
    1582: .    vu - input vector (ghosted)
    1583: .    J - output matrix
    1584: -    w - any user data

    1586:     Notes: Does NOT do ghost updates on vu upon entry

    1588:     Level: advanced

    1590: .seealso: DAFormFunction1()

    1592: @*/
    1593: PetscErrorCode PETSCDM_DLLEXPORT DAComputeJacobian1(DA da,Vec vu,Mat J,void *w)
    1594: {
    1596:   void           *u;
    1597:   DALocalInfo    info;

    1600:   DAGetLocalInfo(da,&info);
    1601:   DAVecGetArray(da,vu,&u);
    1602:   (*da->lj)(&info,u,J,w);
    1603:   DAVecRestoreArray(da,vu,&u);
    1604:   return(0);
    1605: }

    1610: /*
    1611:     DAComputeJacobian1WithAdifor - Evaluates a ADIFOR provided Jacobian local function on each processor that
    1612:         share a DA

    1614:    Input Parameters:
    1615: +    da - the DA that defines the grid
    1616: .    vu - input vector (ghosted)
    1617: .    J - output matrix
    1618: -    w - any user data

    1620:     Notes: Does NOT do ghost updates on vu upon entry

    1622: .seealso: DAFormFunction1()

    1624: */
    1625: PetscErrorCode PETSCDM_DLLEXPORT DAComputeJacobian1WithAdifor(DA da,Vec vu,Mat J,void *w)
    1626: {
    1627:   PetscErrorCode  ierr;
    1628:   PetscInt        i,Nc,N;
    1629:   ISColoringValue *color;
    1630:   DALocalInfo     info;
    1631:   PetscScalar     *u,*g_u,*g_f,*f,*p_u;
    1632:   ISColoring      iscoloring;
    1633:   void            (*lf)(PetscInt*,DALocalInfo*,PetscScalar*,PetscScalar*,PetscInt*,PetscScalar*,PetscScalar*,PetscInt*,void*,PetscErrorCode*) =
    1634:                   (void (*)(PetscInt*,DALocalInfo*,PetscScalar*,PetscScalar*,PetscInt*,PetscScalar*,PetscScalar*,PetscInt*,void*,PetscErrorCode*))*da->adifor_lf;

    1637:   DAGetColoring(da,IS_COLORING_GHOSTED,&iscoloring);
    1638:   Nc   = iscoloring->n;
    1639:   DAGetLocalInfo(da,&info);
    1640:   N    = info.gxm*info.gym*info.gzm*info.dof;

    1642:   /* get space for derivative objects.  */
    1643:   PetscMalloc(Nc*info.gxm*info.gym*info.gzm*info.dof*sizeof(PetscScalar),&g_u);
    1644:   PetscMemzero(g_u,Nc*info.gxm*info.gym*info.gzm*info.dof*sizeof(PetscScalar));
    1645:   p_u   = g_u;
    1646:   color = iscoloring->colors;
    1647:   for (i=0; i<N; i++) {
    1648:     p_u[*color++] = 1.0;
    1649:     p_u          += Nc;
    1650:   }
    1651:   ISColoringDestroy(iscoloring);
    1652:   PetscMalloc(Nc*info.xm*info.ym*info.zm*info.dof*sizeof(PetscScalar),&g_f);
    1653:   PetscMalloc(info.xm*info.ym*info.zm*info.dof*sizeof(PetscScalar),&f);

    1655:   /* Seed the input array g_u with coloring information */
    1656:
    1657:   VecGetArray(vu,&u);
    1658:   (lf)(&Nc,&info,u,g_u,&Nc,f,g_f,&Nc,w,&ierr);
    1659:   VecRestoreArray(vu,&u);

    1661:   /* stick the values into the matrix */
    1662:   /* PetscScalarView(Nc*info.xm*info.ym,g_f,0); */
    1663:   MatSetValuesAdifor(J,Nc,g_f);

    1665:   /* return space for derivative objects.  */
    1666:   PetscFree(g_u);
    1667:   PetscFree(g_f);
    1668:   PetscFree(f);
    1669:   return(0);
    1670: }

    1674: /*@C
    1675:     DAMultiplyByJacobian1WithAD - Applies a Jacobian function supplied by ADIFOR or ADIC
    1676:     to a vector on each processor that shares a DA.

    1678:    Input Parameters:
    1679: +    da - the DA that defines the grid
    1680: .    vu - Jacobian is computed at this point (ghosted)
    1681: .    v - product is done on this vector (ghosted)
    1682: .    fu - output vector = J(vu)*v (not ghosted)
    1683: -    w - any user data

    1685:     Notes:
    1686:     This routine does NOT do ghost updates on vu and v upon entry.
    1687:
    1688:     Automatically calls DAMultiplyByJacobian1WithAdifor() or DAMultiplyByJacobian1WithAdic()
    1689:     depending on whether DASetLocalAdicMFFunction() or DASetLocalAdiforMFFunction() was called.

    1691:    Level: advanced

    1693: .seealso: DAFormFunction1(), DAMultiplyByJacobian1WithAdifor(), DAMultiplyByJacobian1WithAdic()

    1695: @*/
    1696: PetscErrorCode PETSCDM_DLLEXPORT DAMultiplyByJacobian1WithAD(DA da,Vec u,Vec v,Vec f,void *w)
    1697: {

    1701:   if (da->adicmf_lf) {
    1702: #if defined(PETSC_HAVE_ADIC)
    1703:     DAMultiplyByJacobian1WithAdic(da,u,v,f,w);
    1704: #else
    1705:     SETERRQ(PETSC_ERR_SUP_SYS,"Requires ADIC to be installed and cannot use complex numbers");
    1706: #endif
    1707:   } else if (da->adiformf_lf) {
    1708:     DAMultiplyByJacobian1WithAdifor(da,u,v,f,w);
    1709:   } else {
    1710:     SETERRQ(PETSC_ERR_ORDER,"Must call DASetLocalAdiforMFFunction() or DASetLocalAdicMFFunction() before using");
    1711:   }
    1712:   return(0);
    1713: }

    1718: /*@C
    1719:     DAMultiplyByJacobian1WithAdifor - Applies a ADIFOR provided Jacobian function on each processor that
    1720:         share a DA to a vector

    1722:    Input Parameters:
    1723: +    da - the DA that defines the grid
    1724: .    vu - Jacobian is computed at this point (ghosted)
    1725: .    v - product is done on this vector (ghosted)
    1726: .    fu - output vector = J(vu)*v (not ghosted)
    1727: -    w - any user data

    1729:     Notes: Does NOT do ghost updates on vu and v upon entry

    1731:    Level: advanced

    1733: .seealso: DAFormFunction1()

    1735: @*/
    1736: PetscErrorCode PETSCDM_DLLEXPORT DAMultiplyByJacobian1WithAdifor(DA da,Vec u,Vec v,Vec f,void *w)
    1737: {
    1739:   PetscScalar    *au,*av,*af,*awork;
    1740:   Vec            work;
    1741:   DALocalInfo    info;
    1742:   void           (*lf)(DALocalInfo*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,void*,PetscErrorCode*) =
    1743:                  (void (*)(DALocalInfo*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,void*,PetscErrorCode*))*da->adiformf_lf;

    1746:   DAGetLocalInfo(da,&info);

    1748:   DAGetGlobalVector(da,&work);
    1749:   VecGetArray(u,&au);
    1750:   VecGetArray(v,&av);
    1751:   VecGetArray(f,&af);
    1752:   VecGetArray(work,&awork);
    1753:   (lf)(&info,au,av,awork,af,w,&ierr);
    1754:   VecRestoreArray(u,&au);
    1755:   VecRestoreArray(v,&av);
    1756:   VecRestoreArray(f,&af);
    1757:   VecRestoreArray(work,&awork);
    1758:   DARestoreGlobalVector(da,&work);

    1760:   return(0);
    1761: }

    1765: /*@C
    1766:        DASetInterpolationType - Sets the type of interpolation that will be
    1767:           returned by DAGetInterpolation()

    1769:    Collective on DA

    1771:    Input Parameter:
    1772: +  da - initial distributed array
    1773: .  ctype - DA_Q1 and DA_Q0 are currently the only supported forms

    1775:    Level: intermediate

    1777:    Notes: you should call this on the coarser of the two DAs you pass to DAGetInterpolation()

    1779: .keywords:  distributed array, interpolation

    1781: .seealso: DACreate1d(), DACreate2d(), DACreate3d(), DADestroy(), DA, DAInterpolationType
    1782: @*/
    1783: PetscErrorCode PETSCDM_DLLEXPORT DASetInterpolationType(DA da,DAInterpolationType ctype)
    1784: {
    1787:   da->interptype = ctype;
    1788:   return(0);
    1789: }