Actual source code: aodatabasic.c

      1: #define PETSCDM_DLL

      3: /*
      4:   The most basic AOData routines. These store the entire database on each processor.
      5:   These routines are very simple; note that we do not even use a private data structure
      6:   for AOData, and the private datastructure for AODataSegment is just used as a simple array.

      8:   These are made slightly complicated by having to be able to handle logical variables
      9:   stored in bit arrays. Thus,
     10:     - Before mallocing to hold a bit array, we shrunk the array length by a factor
     11:       of 8 using PetscBTLength()
     12:     - We use PetscBitMemcpy() to allow us to copy at the individual bit level;
     13:       for regular datatypes this just does a regular memcpy().
     14: */

     16:  #include src/dm/ao/aoimpl.h
     17:  #include petscsys.h
     18:  #include petscbt.h

     22: PetscErrorCode AODataDestroy_Basic(AOData ao)
     23: {
     25:   AODataKey     *key = ao->keys,*nextkey;
     26:   AODataSegment *seg,*nextseg;

     29:   while (key) {
     30:     PetscFree(key->name);
     31:     if (key->ltog) {
     32:       ISLocalToGlobalMappingDestroy(key->ltog);
     33:     }
     34:     seg = key->segments;
     35:     while (seg) {
     36:       PetscFree(seg->data);
     37:       PetscFree(seg->name);
     38:       nextseg = seg->next;
     39:       PetscFree(seg);
     40:       seg     = nextseg;
     41:     }
     42:     PetscFree(key->rowners);
     43:     nextkey = key->next;
     44:     PetscFree(key);
     45:     key     = nextkey;
     46:   }
     47:   PetscHeaderDestroy(ao);
     48:   return(0);
     49: }

     53: PetscErrorCode AODataView_Basic_Binary(AOData ao,PetscViewer viewer)
     54: {
     56:   PetscInt       N;
     57:   int            fd;
     58:   AODataSegment  *segment;
     59:   AODataKey      *key = ao->keys;
     60:   char           paddedname[256];

     64:   PetscViewerBinaryGetDescriptor(viewer,&fd);

     66:   /* write out number of keys */
     67:   PetscBinaryWrite(fd,&ao->nkeys,1,PETSC_INT,PETSC_FALSE);

     69:   while (key) {
     70:     N   = key->N;
     71:     /*
     72:        Write out name of key - use a fixed length for the name in the binary
     73:        file to make seeking easier
     74:     */
     75:     PetscMemzero(paddedname,256*sizeof(char));
     76:     PetscStrncpy(paddedname,key->name,255);
     77:     PetscBinaryWrite(fd,paddedname,256,PETSC_CHAR,PETSC_FALSE);
     78:     /* write out the number of indices */
     79:     PetscBinaryWrite(fd,&key->N,1,PETSC_INT,PETSC_FALSE);
     80:     /* write out number of segments */
     81:     PetscBinaryWrite(fd,&key->nsegments,1,PETSC_INT,PETSC_FALSE);
     82:
     83:     /* loop over segments writing them out */
     84:     segment = key->segments;
     85:     while (segment) {
     86:       /*
     87:          Write out name of segment - use a fixed length for the name in the binary
     88:          file to make seeking easier
     89:       */
     90:       PetscMemzero(paddedname,256*sizeof(char));
     91:       PetscStrncpy(paddedname,segment->name,255);
     92:       PetscBinaryWrite(fd,paddedname,256,PETSC_CHAR,PETSC_FALSE);
     93:       PetscBinaryWrite(fd,&segment->bs,1,PETSC_INT,PETSC_FALSE);
     94:       PetscBinaryWrite(fd,&segment->datatype,1,PETSC_INT,PETSC_FALSE);
     95:       /* write out the data */
     96:       PetscBinaryWrite(fd,segment->data,N*segment->bs,segment->datatype,PETSC_FALSE);
     97:       segment = segment->next;
     98:     }
     99:     key = key->next;
    100:   }

    102:   return(0);
    103: }

    105: /*
    106:       All processors have the same data so processor 0 prints it
    107: */
    110: PetscErrorCode AODataView_Basic_ASCII(AOData ao,PetscViewer viewer)
    111: {
    112:   PetscErrorCode    ierr;
    113:   PetscMPIInt       rank,size;
    114:   PetscInt          j,k,l,nkeys,nsegs,i,N,bs,zero = 0;
    115:   char              **keynames,**segnames,*segvalue;
    116:   const char        *stype,*dt;
    117:   AODataSegment     *segment;
    118:   AODataKey         *key = ao->keys;
    119:   PetscDataType     dtype;
    120:   PetscViewerFormat format;

    123:   MPI_Comm_rank(ao->comm,&rank);
    124:   if (rank) return(0);
    125:   MPI_Comm_size(ao->comm,&size);

    127:   PetscViewerGetFormat(viewer,&format);
    128:   if (format == PETSC_VIEWER_ASCII_INFO) {
    129:     AODataGetInfo(ao,&nkeys,&keynames);
    130:     for (i=0; i<nkeys; i++) {
    131:       AODataKeyGetInfo(ao,keynames[i],&N,0,&nsegs,&segnames);
    132:       PetscViewerASCIIPrintf(viewer,"  %s: (%D)\n",keynames[i],N);
    133:       for (j=0; j<nsegs; j++) {
    134:         AODataSegmentGetInfo(ao,keynames[i],segnames[j],&bs,&dtype);
    135:         stype = PetscDataTypes[dtype];
    136:         if (dtype == PETSC_CHAR) {
    137:           AODataSegmentGet(ao,keynames[i],segnames[j],1,&zero,(void **)&segvalue);
    138:           PetscViewerASCIIPrintf(viewer,"      %s: (%D) %s -> %s\n",segnames[j],bs,stype,segvalue);
    139:           AODataSegmentRestore(ao,keynames[i],segnames[j],1,&zero,(void **)&segvalue);
    140:         } else {
    141:           PetscViewerASCIIPrintf(viewer,"      %s: (%D) %s\n",segnames[j],bs,stype);
    142:         }
    143:       }
    144:     }
    145:     PetscFree(keynames);
    146:   } else {
    147:     PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
    148:     while (key) {
    149:       PetscViewerASCIIPrintf(viewer,"AOData Key: %s Length %D Ownership: ",key->name,key->N);
    150:       for (j=0; j<size+1; j++) {PetscViewerASCIIPrintf(viewer,"%D ",key->rowners[j]);}
    151:       PetscViewerASCIIPrintf(viewer,"\n");

    153:       segment = key->segments;
    154:       while (segment) {
    155:         dt = PetscDataTypes[segment->datatype];
    156:         PetscViewerASCIIPrintf(viewer,"  AOData Segment: %s Blocksize %D datatype %s\n",segment->name,segment->bs,dt);
    157:         if (segment->datatype == PETSC_INT) {
    158:           PetscInt *mdata = (PetscInt*)segment->data;
    159:           for (k=0; k<key->N; k++) {
    160:             PetscViewerASCIIPrintf(viewer," %D: ",k);
    161:             for (l=0; l<segment->bs; l++) {
    162:               PetscViewerASCIIPrintf(viewer,"   %D ",mdata[k*segment->bs + l]);
    163:             }
    164:             PetscViewerASCIIPrintf(viewer,"\n");
    165:           }
    166:         } else if (segment->datatype == PETSC_DOUBLE) {
    167:           PetscReal *mdata = (PetscReal*)segment->data;
    168:           for (k=0; k<key->N; k++) {
    169:             PetscViewerASCIIPrintf(viewer," %D: ",k);
    170:             for (l=0; l<segment->bs; l++) {
    171:               PetscViewerASCIIPrintf(viewer,"   %18.16e ",mdata[k*segment->bs + l]);
    172:             }
    173:             PetscViewerASCIIPrintf(viewer,"\n");
    174:           }
    175:         } else if (segment->datatype == PETSC_SCALAR) {
    176:           PetscScalar *mdata = (PetscScalar*)segment->data;
    177:           for (k=0; k<key->N; k++) {
    178:             PetscViewerASCIIPrintf(viewer," %D: ",k);
    179:             for (l=0; l<segment->bs; l++) {
    180: #if !defined(PETSC_USE_COMPLEX)
    181:               PetscViewerASCIIPrintf(viewer,"   %18.16e ",mdata[k*segment->bs + l]);
    182: #else
    183:               PetscScalar x = mdata[k*segment->bs + l];
    184:               if (PetscImaginaryPart(x) > 0.0) {
    185:                 PetscViewerASCIIPrintf(viewer," %18.16e + %18.16e i \n",PetscRealPart(x),PetscImaginaryPart(x));
    186:               } else if (PetscImaginaryPart(x) < 0.0) {
    187:                 PetscViewerASCIIPrintf(viewer,"   %18.16e - %18.16e i \n",PetscRealPart(x),-PetscImaginaryPart(x));
    188:               } else {
    189:                 PetscViewerASCIIPrintf(viewer,"   %18.16e \n",PetscRealPart(x));
    190:               }
    191: #endif
    192:             }
    193:           }
    194:           PetscViewerASCIIPrintf(viewer,"\n");
    195:         } else if (segment->datatype == PETSC_LOGICAL) {
    196:           PetscBT mdata = (PetscBT) segment->data;
    197:           for (k=0; k<key->N; k++) {
    198:             PetscViewerASCIIPrintf(viewer," %D: ",k);
    199:             for (l=0; l<segment->bs; l++) {
    200:               PetscViewerASCIIPrintf(viewer,"   %d ",(int)PetscBTLookup(mdata,k*segment->bs + l));
    201:             }
    202:             PetscViewerASCIIPrintf(viewer,"\n");
    203:           }
    204:         } else if (segment->datatype == PETSC_CHAR) {
    205:           char * mdata = (char*)segment->data;
    206:           for (k=0; k<key->N; k++) {
    207:             PetscViewerASCIIPrintf(viewer,"  %s ",mdata + k*segment->bs);
    208:           }
    209:           PetscViewerASCIIPrintf(viewer,"\n");
    210:         } else {
    211:           SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Unknown PETSc data format");
    212:         }
    213:         segment = segment->next;
    214:       }
    215:       key = key->next;
    216:     }
    217:   }
    218:   PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
    219:   return(0);
    220: }

    224: PetscErrorCode AODataView_Basic(AOData ao,PetscViewer viewer)
    225: {
    227:   PetscMPIInt    rank;
    228:   PetscTruth     iascii,isbinary;

    231:   MPI_Comm_rank(ao->comm,&rank);
    232:   if (rank) return(0);

    234:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    235:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);
    236:   if (iascii) {
    237:     AODataView_Basic_ASCII(ao,viewer);
    238:   } else if (isbinary) {
    239:     AODataView_Basic_Binary(ao,viewer);
    240:   } else {
    241:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for AOData basic",((PetscObject)viewer)->type_name);
    242:   }

    244:   return(0);
    245: }

    249: PetscErrorCode AODataKeyRemove_Basic(AOData aodata,const char name[])
    250: {
    251:   AODataSegment    *segment,*iseg;
    252:   AODataKey        *key,*ikey;
    254:   PetscTruth       flag;

    257:   AODataKeyFind_Private(aodata,name,&flag,&key);
    258:   if (flag != 1)  return(0);
    259:   aodata->nkeys--;

    261:   segment = key->segments;
    262:   while (segment) {
    263:     iseg    = segment->next;
    264:     PetscFree(segment->name);
    265:     PetscFree(segment->data);
    266:     PetscFree(segment);
    267:     segment = iseg;
    268:   }
    269:   ikey = aodata->keys;
    270:   if (key == ikey) {
    271:     aodata->keys = key->next;
    272:     goto finishup1;
    273:   }
    274:   while (ikey->next != key) {
    275:     ikey = ikey->next;
    276:   }
    277:   ikey->next = key->next;

    279:   finishup1:

    281:   PetscFree(key->name);
    282:   PetscFree(key->rowners);
    283:   PetscFree(key);

    285:   return(0);
    286: }

    290: PetscErrorCode AODataSegmentRemove_Basic(AOData aodata,const char name[],const char segname[])
    291: {
    292:   AODataSegment    *segment,*iseg;
    293:   AODataKey        *key;
    295:   PetscTruth       flag;

    298:   AODataSegmentFind_Private(aodata,name,segname,&flag,&key,&iseg);
    299:   if (!flag)  return(0);
    300:   key->nsegments--;

    302:   segment = key->segments;
    303:   if (segment == iseg) {
    304:     key->segments = segment->next;
    305:     goto finishup2;
    306:   }
    307:   while (segment->next != iseg) {
    308:     segment = segment->next;
    309:   }
    310:   segment->next = iseg->next;
    311:   segment       = iseg;
    312:
    313:   finishup2:

    315:   PetscFree(segment->name);
    316:   PetscFree(segment->data);
    317:   PetscFree(segment);
    318:   return(0);
    319: }

    324: PetscErrorCode AODataSegmentAdd_Basic(AOData aodata,const char name[],const char segname[],PetscInt bs,PetscInt n,PetscInt *keys,void *data,PetscDataType dtype)
    325: {
    326:   AODataSegment  *segment,*iseg;
    327:   AODataKey      *key;
    329:   PetscInt       N,i,*akeys,datasize,*fkeys,Nb,j;
    330:   PetscMPIInt    size,*lens,*disp,nn;
    331:   MPI_Datatype   mtype;
    332:   char           *adata;
    333:   MPI_Comm       comm = aodata->comm;
    334:   PetscTruth     flag;

    337:   AODataKeyFind_Private(aodata,name,&flag,&key);
    338:   if (!flag)  SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Key %s doesn't exist",name);
    339:   AODataSegmentFind_Private(aodata,name,segname,&flag,&key,&iseg);
    340:   if (flag)  SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Segment %s in key %s already exists",name,segname);

    342:   PetscNew(AODataSegment,&segment);
    343:   if (iseg) {
    344:     iseg->next    = segment;
    345:   } else {
    346:     key->segments = segment;
    347:   }
    348:   segment->next     = 0;
    349:   segment->bs       = bs;
    350:   segment->datatype = dtype;

    352:   PetscDataTypeGetSize(dtype,&datasize);

    354:   /*
    355:      If keys not given, assume each processor provides entire data
    356:   */
    357:   if (!keys && n == key->N) {
    358:     char *fdata1;
    359:     if (dtype == PETSC_LOGICAL) Nb = PetscBTLength(key->N); else Nb = key->N;
    360:     PetscMalloc((Nb*bs+1)*datasize,&fdata1);
    361:     PetscBitMemcpy(fdata1,0,data,0,key->N*bs,dtype);
    362:     segment->data = (void*)fdata1;
    363:   } else if (!keys) {
    364:     SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Keys not given, but not all data given on each processor");
    365:   } else {
    366:     /* transmit all lengths to all processors */
    367:     MPI_Comm_size(comm,&size);
    368:     PetscMalloc(2*size*sizeof(PetscMPIInt),&lens);
    369:     disp = lens + size;
    370:     nn   = n;
    371:     MPI_Allgather(&nn,1,MPI_INT,lens,1,MPI_INT,comm);
    372:     N =  0;
    373:     for (i=0; i<size; i++) {
    374:       disp[i]  = N;
    375:       N       += lens[i];
    376:     }
    377:     if (N != key->N) {
    378:       SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Did not provide correct number of keys for keyname");
    379:     }

    381:     /*
    382:       Allocate space for all keys and all data
    383:     */
    384:     PetscMalloc((N+1)*sizeof(PetscInt),&akeys);
    385:     PetscMalloc((N*bs+1)*datasize,&adata);

    387:     MPI_Allgatherv(keys,n,MPIU_INT,akeys,lens,disp,MPIU_INT,comm);
    388:     for (i=0; i<size; i++) {
    389:       disp[i] *= bs;
    390:       lens[i] *= bs;
    391:     }
    392:
    393:     if (dtype != PETSC_LOGICAL) {
    394:       char *fdata2;

    396:       PetscDataTypeToMPIDataType(dtype,&mtype);
    397:       MPI_Allgatherv(data,n*bs,mtype,adata,lens,disp,mtype,comm);
    398:       PetscFree(lens);

    400:       /*
    401:         Now we have all the keys and data we need to put it in order
    402:       */
    403:       PetscMalloc((key->N+1)*sizeof(PetscInt),&fkeys);
    404:       PetscMemzero(fkeys,(key->N+1)*sizeof(PetscInt));
    405:       PetscMalloc((key->N*bs+1)*datasize,&fdata2);

    407:       for (i=0; i<N; i++) {
    408:         if (fkeys[akeys[i]] != 0) {
    409:           SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Duplicate key");
    410:         }
    411:         if (fkeys[akeys[i]] >= N) {
    412:           SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Key out of range");
    413:         }
    414:         fkeys[akeys[i]] = 1;
    415:         PetscBitMemcpy(fdata2,akeys[i]*bs,adata,i*bs,bs,dtype);
    416:       }
    417:       for (i=0; i<N; i++) {
    418:         if (!fkeys[i]) {
    419:           SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Missing key");
    420:         }
    421:       }
    422:       segment->data = (void*)fdata2;
    423:     } else {
    424:       /*
    425:             For logical input the length is given by the user in bits; we need to
    426:             convert to bytes to send with MPI
    427:       */
    428:       PetscBT fdata3;
    429:       PetscBT mvalues = (PetscBT) data;
    430:       char *values;
    431:       PetscMalloc((n+1)*bs*sizeof(char),&values);
    432:       for (i=0; i<n; i++) {
    433:         for (j=0; j<bs; j++) {
    434:           if (PetscBTLookup(mvalues,i*bs+j)) values[i*bs+j] = 1; else values[i*bs+j] = 0;
    435:         }
    436:       }

    438:       MPI_Allgatherv(values,n*bs,MPI_BYTE,adata,lens,disp,MPI_BYTE,comm);
    439:       PetscFree(lens);
    440:       PetscFree(values);

    442:       /*
    443:         Now we have all the keys and data we need to put it in order
    444:       */
    445:       PetscMalloc((key->N+1)*sizeof(PetscInt),&fkeys);
    446:       PetscMemzero(fkeys,(key->N+1)*sizeof(PetscInt));
    447:       PetscBTCreate(N*bs,fdata3);

    449:       for (i=0; i<N; i++) {
    450:         if (fkeys[akeys[i]] != 0) {
    451:           SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Duplicate key");
    452:         }
    453:         if (fkeys[akeys[i]] >= N) {
    454:           SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Key out of range");
    455:         }
    456:         fkeys[akeys[i]] = 1;
    457:         for (j=0; j<bs; j++) {
    458:           if (adata[i*bs+j]) { PetscBTSet(fdata3,i*bs+j); }
    459:         }
    460:       }
    461:       for (i=0; i<N; i++) {
    462:         if (!fkeys[i]) {
    463:           SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Missing key");
    464:         }
    465:       }
    466:       segment->data = (void*)fdata3;
    467:     }
    468:     PetscFree(akeys);
    469:     PetscFree(adata);
    470:     PetscFree(fkeys);
    471:   }

    473:   key->nsegments++;

    475:   PetscStrallocpy(segname,&segment->name);
    476:   return(0);
    477: }

    481: PetscErrorCode AODataSegmentGetExtrema_Basic(AOData ao,const char name[],const char segname[],void *xmax,void *xmin)
    482: {
    483:   AODataSegment  *segment;
    484:   AODataKey      *key;
    486:   PetscInt       i,bs,n,j;
    487:   PetscTruth     flag;

    490:   /* find the correct segment */
    491:   AODataSegmentFind_Private(ao,name,segname,&flag,&key,&segment);
    492:   if (!flag) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Cannot locate segment");

    494:   n       = key->N;
    495:   bs      = segment->bs;

    497:   if (segment->datatype == PETSC_INT) {
    498:     PetscInt *vmax = (PetscInt*)xmax,*vmin = (PetscInt*)xmin,*values = (PetscInt*)segment->data;
    499:     for (j=0; j<bs; j++) {
    500:       vmax[j] = vmin[j] = values[j];
    501:     }
    502:     for (i=1; i<n; i++) {
    503:       for (j=0; j<bs; j++) {
    504:         vmax[j] = PetscMax(vmax[j],values[bs*i+j]);
    505:         vmin[j] = PetscMin(vmin[j],values[bs*i+j]);
    506:       }
    507:     }
    508:   } else if (segment->datatype == PETSC_DOUBLE) {
    509:     PetscReal *vmax = (PetscReal*)xmax,*vmin = (PetscReal*)xmin,*values = (PetscReal*)segment->data;
    510:     for (j=0; j<bs; j++) {
    511:       vmax[j] = vmin[j] = values[j];
    512:     }
    513:     for (i=1; i<n; i++) {
    514:       for (j=0; j<bs; j++) {
    515:         vmax[j] = PetscMax(vmax[j],values[bs*i+j]);
    516:         vmin[j] = PetscMin(vmin[j],values[bs*i+j]);
    517:       }
    518:     }
    519:   } else SETERRQ(PETSC_ERR_SUP,"Cannot find extrema for this data type");

    521:   return(0);
    522: }

    526: PetscErrorCode AODataSegmentGet_Basic(AOData ao,const char name[],const char segname[],PetscInt n,PetscInt *keys,void **data)
    527: {
    528:   AODataSegment  *segment;
    529:   AODataKey      *key;
    531:   PetscInt       dsize,i,bs,nb;
    532:   char           *idata,*odata;
    533:   PetscTruth     flag;
    534:
    536:   /* find the correct segment */
    537:   AODataSegmentFind_Private(ao,name,segname,&flag,&key,&segment);
    538:   if (!flag) SETERRQ2(PETSC_ERR_ARG_WRONG,"Cannot locate segment %s in key %s",segname,name);

    540:   PetscDataTypeGetSize(segment->datatype,&dsize);
    541:   bs    = segment->bs;
    542:   if (segment->datatype == PETSC_LOGICAL) nb = PetscBTLength(n); else nb = n;
    543:   PetscMalloc((nb+1)*bs*dsize,&odata);
    544:   idata = (char*)segment->data;
    545:   for (i=0; i<n; i++) {
    546:     PetscBitMemcpy(odata,i*bs,idata,keys[i]*bs,bs,segment->datatype);
    547:   }
    548:   *data = (void*)odata;
    549:   return(0);
    550: }

    554: PetscErrorCode AODataSegmentRestore_Basic(AOData aodata,const char name[],const char segname[],PetscInt n,PetscInt *keys,void **data)
    555: {

    559:   PetscFree(*data);
    560:   return(0);
    561: }

    565: PetscErrorCode AODataSegmentGetLocal_Basic(AOData ao,const char name[],const char segname[],PetscInt n,PetscInt *keys,void **data)
    566: {
    568:   PetscInt       *globals,*locals,bs;
    569:   PetscDataType  dtype;
    570:   AODataKey      *key;
    571:   PetscTruth     flag;

    574:   AODataKeyFind_Private(ao,segname,&flag,&key);
    575:   if (!flag) SETERRQ(PETSC_ERR_ARG_WRONG,"Segment does not have corresponding key");
    576:   if (!key->ltog) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"No local to global mapping set for key");
    577:   AODataSegmentGetInfo(ao,name,segname,&bs,&dtype);
    578:   if (dtype != PETSC_INT) SETERRQ(PETSC_ERR_ARG_WRONG,"Datatype of segment must be PETSC_INT");

    580:   /* get the values in global indexing */
    581:   AODataSegmentGet_Basic(ao,name,segname,n,keys,(void **)&globals);
    582:
    583:   /* allocate space to store them in local indexing */
    584:   PetscMalloc((n+1)*bs*sizeof(PetscInt),&locals);

    586:   ISGlobalToLocalMappingApply(key->ltog,IS_GTOLM_MASK,n*bs,globals,PETSC_NULL,locals);

    588:   AODataSegmentRestore_Basic(ao,name,segname,n,keys,(void **)&globals);

    590:   *data = (void*)locals;
    591:   return(0);
    592: }

    596: PetscErrorCode AODataSegmentRestoreLocal_Basic(AOData aodata,const char name[],const char segname[],PetscInt n,PetscInt *keys,void **data)
    597: {

    601:   PetscFree(*data);
    602:   return(0);
    603: }

    605: EXTERN PetscErrorCode AOBasicGetIndices_Private(AO,PetscInt **,PetscInt **);

    609: PetscErrorCode AODataKeyRemap_Basic(AOData aodata,const char keyname[],AO ao)
    610: {
    612:   PetscInt       *inew,k,*ii,nk,dsize,bs,nkb;
    613:   char           *data,*tmpdata;
    614:   AODataKey      *key;
    615:   AODataSegment  *seg;
    616:   PetscTruth     flag,match;

    620:   /* remap all the values in the segments that match the key */
    621:   key = aodata->keys;
    622:   while (key) {
    623:     seg = key->segments;
    624:     while (seg) {
    625:       PetscStrcmp(seg->name,keyname,&match);
    626:       if (!match) {
    627:         seg = seg->next;
    628:         continue;
    629:       }
    630:       if (seg->datatype != PETSC_INT) {
    631:         SETERRQ(PETSC_ERR_ARG_WRONG,"Segment name same as key but not integer type");
    632:       }
    633:       nk   = seg->bs*key->N;
    634:       ii   = (PetscInt*)seg->data;
    635:       AOPetscToApplication(ao,nk,ii);
    636:       seg  = seg->next;
    637:     }
    638:     key = key->next;
    639:   }
    640:
    641:   AOBasicGetIndices_Private(ao,&inew,PETSC_NULL);
    642:   /* reorder in the arrays all the values for the key */
    643:   AODataKeyFind_Private(aodata,keyname,&flag,&key);
    644:   if (!flag) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Could not find key");
    645:   nk  = key->N;
    646:   seg = key->segments;
    647:   while (seg) {
    648:     PetscDataTypeGetSize(seg->datatype,&dsize);
    649:     bs      = seg->bs;
    650:     data    = (char*)seg->data;
    651:     if (seg->datatype == PETSC_LOGICAL) nkb = PetscBTLength(nk*bs); else nkb = nk*bs;
    652:     PetscMalloc((nkb+1)*dsize,&tmpdata);

    654:     for (k=0; k<nk; k++) {
    655:       PetscBitMemcpy(tmpdata,inew[k]*bs,data,k*bs,bs,seg->datatype);
    656:     }
    657:     PetscMemcpy(data,tmpdata,nkb*dsize);
    658:     PetscFree(tmpdata);
    659:     seg = seg->next;
    660:   }

    662:   return(0);
    663: }

    667: PetscErrorCode AODataKeyGetAdjacency_Basic(AOData aodata,const char keyname[],Mat *adj)
    668: {
    670:   PetscInt       cnt,i,j,*jj,*ii,nlocal,n,*nb,bs,ls;
    671:   AODataKey     *key;
    672:   AODataSegment *seg;
    673:   PetscTruth    flag;

    676:   AODataSegmentFind_Private(aodata,keyname,keyname,&flag,&key,&seg);
    677:   if (!flag) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Cannot locate key with neighbor segment");

    679:   /*
    680:      Get the beginning of the neighbor list for this processor
    681:   */
    682:   bs     = seg->bs;
    683:   nb     = (PetscInt*)seg->data;
    684:   nb    += bs*key->rstart;
    685:   nlocal = key->rend - key->rstart;
    686:   n      = bs*key->N;

    688:   /*
    689:       Assemble the adjacency graph: first we determine total number of entries
    690:   */
    691:   cnt = 0;
    692:   for (i=0; i<bs*nlocal; i++) {
    693:     if (nb[i] >= 0) cnt++;
    694:   }
    695:   PetscMalloc((nlocal + 1)*sizeof(PetscInt),&ii);
    696:   PetscMalloc((cnt+1)*sizeof(PetscInt),&jj);
    697:   ii[0] = 0;
    698:   cnt   = 0;
    699:   for (i=0; i<nlocal; i++) {
    700:     ls = 0;
    701:     for (j=0; j<bs; j++) {
    702:       if (nb[bs*i+j] >= 0) {
    703:         jj[cnt++] = nb[bs*i+j];
    704:         ls++;
    705:       }
    706:     }
    707:     /* now sort the column indices for this row */
    708:     PetscSortInt(ls,jj+cnt-ls);
    709:     ii[i+1] = cnt;
    710:   }

    712:   MatCreate(aodata->comm,adj);
    713:   MatSetSizes(*adj,nlocal,n,PETSC_DETERMINE,n);
    714:   MatSetType(*adj,MATMPIADJ);
    715:   MatMPIAdjSetPreallocation(*adj,ii,jj,PETSC_NULL);
    716:   return(0);
    717: }

    721: PetscErrorCode AODataSegmentPartition_Basic(AOData aodata,const char keyname[],const char segname[])
    722: {
    724:   PetscMPIInt    size;
    725:   PetscInt        bs,i,j,*idx,nc,*isc;
    726:   AO             ao;
    727:   AODataKey      *key,*keyseg;
    728:   AODataSegment  *segment;
    729:   PetscTruth     flag;

    733:   AODataKeyFind_Private(aodata,segname,&flag,&keyseg);
    734:   if (!flag) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Cannot locate segment as a key");
    735:   PetscMalloc(keyseg->N*sizeof(PetscInt),&isc);
    736:   PetscMemzero(isc,keyseg->N*sizeof(PetscInt));

    738:   AODataSegmentFind_Private(aodata,keyname,segname,&flag,&key,&segment);
    739:   if (flag != 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Cannot locate segment");
    740:   MPI_Comm_size(aodata->comm,&size);

    742:   bs                = segment->bs;

    744:   idx = (PetscInt*)segment->data;
    745:   nc  = 0;
    746:   for (i=0; i<size; i++) {
    747:     for (j=bs*key->rowners[i]; j<bs*key->rowners[i+1]; j++) {
    748:       /* allow some keys to have fewer data than others, indicate with a -1 */
    749:       if (idx[j] >= 0) {
    750:         if (!isc[idx[j]]) {
    751:           isc[idx[j]] = ++nc;
    752:         }
    753:       }
    754:     }
    755:   }
    756:   for (i=0; i<nc; i++) {
    757:     isc[i]--;
    758:   }

    760:   AOCreateBasic(aodata->comm,keyseg->nlocal,isc+keyseg->rstart,PETSC_NULL,&ao);
    761:   PetscFree(isc);

    763:   AODataKeyRemap(aodata,segname,ao);
    764:   AODestroy(ao);
    765:   return(0);
    766: }

    770: PetscErrorCode AODataKeyGetActive_Basic(AOData aodata,const char name[],const char segname[],PetscInt n,PetscInt *keys,PetscInt wl,IS *is)
    771: {
    773:   PetscInt       i,cnt,*fnd,bs;
    774:   AODataKey      *key;
    775:   AODataSegment  *segment;
    776:   PetscBT        bt;
    777:   PetscTruth     flag;

    780:   AODataSegmentFind_Private(aodata,name,segname,&flag,&key,&segment);
    781:   if (!flag) SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot locate segment");

    783:   bt = (PetscBT) segment->data;
    784:   bs = segment->bs;

    786:   if (wl >= bs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Bit field (wl) argument too large");

    788:   /* count the active ones */
    789:   cnt = 0;
    790:   for (i=0; i<n; i++) {
    791:     if (PetscBTLookup(bt,keys[i]*bs+wl)) {
    792:       cnt++;
    793:     }
    794:   }

    796:   PetscMalloc((cnt+1)*sizeof(PetscInt),&fnd);
    797:   cnt = 0;
    798:   for (i=0; i<n; i++) {
    799:     if (PetscBTLookup(bt,keys[i]*bs+wl)) {
    800:       fnd[cnt++] = keys[i];
    801:     }
    802:   }
    803:
    804:   ISCreateGeneral(aodata->comm,cnt,fnd,is);
    805:   PetscFree(fnd);
    806:   return(0);
    807: }

    811: PetscErrorCode AODataKeyGetActiveLocal_Basic(AOData aodata,const char name[],const char segname[],PetscInt n,PetscInt *keys,PetscInt wl,IS *is)
    812: {
    814:   PetscInt       i,cnt,*fnd,bs,*locals;
    815:   AODataKey      *key;
    816:   AODataSegment  *segment;
    817:   PetscBT        bt;
    818:   PetscTruth     flag;

    821:   AODataSegmentFind_Private(aodata,name,segname,&flag,&key,&segment);
    822:   if (!flag) SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot locate segment");

    824:   bt = (PetscBT) segment->data;
    825:   bs = segment->bs;

    827:   if (wl >= bs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Bit field (wl) argument too large");

    829:   /* count the active ones */
    830:   cnt = 0;
    831:   for (i=0; i<n; i++) {
    832:     if (PetscBTLookup(bt,keys[i]*bs+wl)) {
    833:       cnt++;
    834:     }
    835:   }

    837:   PetscMalloc((cnt+1)*sizeof(PetscInt),&fnd);
    838:   cnt = 0;
    839:   for (i=0; i<n; i++) {
    840:     if (PetscBTLookup(bt,keys[i]*bs+wl)) {
    841:       fnd[cnt++] = keys[i];
    842:     }
    843:   }
    844:
    845:   PetscMalloc((n+1)*sizeof(PetscInt),&locals);
    846:   ISGlobalToLocalMappingApply(key->ltog,IS_GTOLM_MASK,cnt,fnd,PETSC_NULL,locals);
    847:   PetscFree(fnd);
    848:   ISCreateGeneral(aodata->comm,cnt,locals,is);
    849:   PetscFree(locals);
    850:   return(0);
    851: }

    853: EXTERN PetscErrorCode AODataSegmentGetReduced_Basic(AOData,const char[],const char[],PetscInt,PetscInt*,IS *);
    854: EXTERN PetscErrorCode AODataPublish_Petsc(PetscObject);

    856: static struct _AODataOps myops = {AODataSegmentAdd_Basic,
    857:                                   AODataSegmentGet_Basic,
    858:                                   AODataSegmentRestore_Basic,
    859:                                   AODataSegmentGetLocal_Basic,
    860:                                   AODataSegmentRestoreLocal_Basic,
    861:                                   AODataSegmentGetReduced_Basic,
    862:                                   AODataSegmentGetExtrema_Basic,
    863:                                   AODataKeyRemap_Basic,
    864:                                   AODataKeyGetAdjacency_Basic,
    865:                                   AODataKeyGetActive_Basic,
    866:                                   AODataKeyGetActiveLocal_Basic,
    867:                                   AODataSegmentPartition_Basic,
    868:                                   AODataKeyRemove_Basic,
    869:                                   AODataSegmentRemove_Basic,
    870:                                   AODataDestroy_Basic,
    871:                                   AODataView_Basic};

    875: /*@C
    876:    AODataCreateBasic - Creates an AO datastructure.

    878:    Collective on MPI_Comm

    880:    Input Parameters:
    881: .  comm  - MPI communicator that is to share AO

    883:    Output Parameter:
    884: .  aoout - the new database

    886:    Options Database Keys:
    887: +  -ao_data_view - Prints entire database at the conclusion of AODataSegmentAdd()
    888: -  -ao_data_view_info - Prints info about database at the conclusion of AODataSegmentAdd()

    890:    Level: intermediate

    892: .keywords: AOData, create

    894: .seealso: AODataSegmentAdd(), AODataDestroy()
    895: @*/
    896: PetscErrorCode PETSCDM_DLLEXPORT AODataCreateBasic(MPI_Comm comm,AOData *aoout)
    897: {
    898:   AOData         ao;

    903:   *aoout = 0;
    904: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
    905:   DMInitializePackage(PETSC_NULL);
    906: #endif

    908:   PetscHeaderCreate(ao,_p_AOData,struct _AODataOps,AODATA_COOKIE,AODATA_BASIC,"AOData",comm,AODataDestroy,AODataView);
    909:   PetscLogObjectMemory(ao,sizeof(struct _p_AOData));

    911:   PetscMemcpy(ao->ops,&myops,sizeof(myops));
    912:   ao->bops->publish = AODataPublish_Petsc;

    914:   ao->nkeys        = 0;
    915:   ao->keys         = 0;
    916:   ao->datacomplete = 0;

    918:   PetscPublishAll(ao);
    919:   *aoout = ao;
    920:   return(0);
    921: }

    925: /*@C
    926:    AODataLoadBasic - Loads an AO database from a file.

    928:    Collective on PetscViewer

    930:    Input Parameters:
    931: .  viewer - the binary file containing the data

    933:    Output Parameter:
    934: .  aoout - the new database

    936:    Options Database Keys:
    937: +  -ao_data_view - Prints entire database at the conclusion of AODataLoadBasic()
    938: -  -ao_data_view_info - Prints info about database at the conclusion of AODataLoadBasic()

    940:    Level: intermediate

    942: .keywords: AOData, create, load, basic

    944: .seealso: AODataSegmentAdd(), AODataDestroy(), AODataCreateBasic(), AODataView()
    945: @*/
    946: PetscErrorCode PETSCDM_DLLEXPORT AODataLoadBasic(PetscViewer viewer,AOData *aoout)
    947: {
    948:   AOData        ao;
    950:   int            fd;
    951:   PetscMPIInt    size,rank;
    952:   PetscInt       nkeys,i,dsize,j,Nb;
    953:   char           paddedname[256];
    954:   AODataSegment  *seg = 0;
    955:   AODataKey      *key = 0;
    956:   MPI_Comm       comm;
    957:   PetscTruth     isbinary,flg1;

    960:   *aoout = 0;
    961: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
    962:   DMInitializePackage(PETSC_NULL);
    963: #endif
    964:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);
    965:   if (!isbinary) {
    966:     SETERRQ(PETSC_ERR_ARG_WRONG,"Viewer must be obtained from PetscViewerBinaryOpen()");
    967:   }

    969:   PetscObjectGetComm((PetscObject)viewer,&comm);
    970:   MPI_Comm_size(comm,&size);
    971:   MPI_Comm_rank(comm,&rank);

    973:   PetscViewerBinaryGetDescriptor(viewer,&fd);

    975:   /* read in number of segments */
    976:   PetscBinaryRead(fd,&nkeys,1,PETSC_INT);

    978:   PetscHeaderCreate(ao,_p_AOData,struct _AODataOps,AODATA_COOKIE,AODATA_BASIC,"AOData",comm,AODataDestroy,AODataView);
    979:   PetscLogObjectMemory(ao,sizeof(struct _p_AOData) + nkeys*sizeof(void*));

    981:   PetscMemcpy(ao->ops,&myops,sizeof(myops));
    982:   ao->bops->publish  = AODataPublish_Petsc;

    984:   ao->nkeys      = nkeys;

    986:   for (i=0; i<nkeys; i++) {
    987:     if (!i) {
    988:       PetscNew(AODataKey,&key);
    989:       ao->keys = key;
    990:     } else {
    991:       PetscNew(AODataKey,&key->next);
    992:       key  = key->next;
    993:     }
    994:     key->ltog = 0;
    995:     key->next = 0;

    997:     /* read in key name */
    998:     PetscBinaryRead(fd,paddedname,256,PETSC_CHAR);
    999:     PetscStrallocpy(paddedname,&key->name);
    1000:     PetscBinaryRead(fd,&key->N,1,PETSC_INT);

    1002:     /* determine Nlocal and rowners for key */
    1003:     key->nlocal  = key->N/size + ((key->N % size) > rank);
    1004:     PetscMalloc((size+1)*sizeof(PetscInt),&key->rowners);
    1005:     MPI_Allgather(&key->nlocal,1,MPI_INT,key->rowners+1,1,MPI_INT,comm);
    1006:     key->rowners[0] = 0;
    1007:     for (j=2; j<=size; j++) {
    1008:       key->rowners[j] += key->rowners[j-1];
    1009:     }
    1010:     key->rstart        = key->rowners[rank];
    1011:     key->rend          = key->rowners[rank+1];

    1013:     /* loop key's segments, reading them in */
    1014:     PetscBinaryRead(fd,&key->nsegments,1,PETSC_INT);

    1016:     for (j=0; j<key->nsegments; j++) {
    1017:       if (!j) {
    1018:         PetscNew(AODataSegment,&seg);
    1019:         key->segments = seg;
    1020:       } else {
    1021:         PetscNew(AODataSegment,&seg->next);
    1022:         seg  = seg->next;
    1023:       }

    1025:       /* read in segment name */
    1026:       PetscBinaryRead(fd,paddedname,256,PETSC_CHAR);
    1027:       PetscStrallocpy(paddedname,&seg->name);

    1029:       /* read in segment blocksize and datatype */
    1030:       PetscBinaryRead(fd,&seg->bs,1,PETSC_INT);
    1031:       PetscBinaryRead(fd,&seg->datatype,1,PETSC_INT);

    1033:       /* allocate the space for the data */
    1034:       PetscDataTypeGetSize(seg->datatype,&dsize);
    1035:       if (seg->datatype == PETSC_LOGICAL) Nb = PetscBTLength(key->N*seg->bs); else Nb = key->N*seg->bs;
    1036:       PetscMalloc(Nb*dsize,&seg->data);
    1037:       /* read in the data */
    1038:       PetscBinaryRead(fd,seg->data,key->N*seg->bs,seg->datatype);
    1039:       seg->next = 0;
    1040:     }
    1041:   }
    1042:   *aoout = ao;

    1044:   PetscOptionsHasName(PETSC_NULL,"-ao_data_view",&flg1);
    1045:   if (flg1) {
    1046:     AODataView(ao,PETSC_VIEWER_STDOUT_(comm));
    1047:   }
    1048:   PetscOptionsHasName(PETSC_NULL,"-ao_data_view_info",&flg1);
    1049:   if (flg1) {
    1050:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(comm),PETSC_VIEWER_ASCII_INFO);
    1051:     AODataView(ao,PETSC_VIEWER_STDOUT_(comm));
    1052:     PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(comm));
    1053:   }
    1054:   PetscPublishAll(ao);
    1055:   return(0);
    1056: }