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: }