Actual source code: general.c

petsc-master 2019-12-15
Report Typos and Errors

  2: /*
  3:      Provides the functions for index sets (IS) defined by a list of integers.
  4: */
  5:  #include <../src/vec/is/is/impls/general/general.h>
  6:  #include <petsc/private/viewerimpl.h>
  7: #include <petsc/private/viewerhdf5impl.h>

  9: static PetscErrorCode ISDuplicate_General(IS is,IS *newIS)
 10: {
 12:   IS_General     *sub = (IS_General*)is->data;
 13:   PetscInt       n;

 16:   PetscLayoutGetLocalSize(is->map, &n);
 17:   ISCreateGeneral(PetscObjectComm((PetscObject) is), n, sub->idx, PETSC_COPY_VALUES, newIS);
 18:   return(0);
 19: }

 21: static PetscErrorCode ISDestroy_General(IS is)
 22: {
 23:   IS_General     *is_general = (IS_General*)is->data;

 27:   if (is_general->allocated) {PetscFree(is_general->idx);}
 28:   PetscObjectComposeFunction((PetscObject)is,"ISGeneralSetIndices_C",NULL);
 29:   PetscObjectComposeFunction((PetscObject)is,"ISGeneralFilter_C",NULL);
 30:   PetscFree(is->data);
 31:   return(0);
 32: }

 34: static PetscErrorCode ISCopy_General(IS is,IS isy)
 35: {
 36:   IS_General     *is_general = (IS_General*)is->data,*isy_general = (IS_General*)isy->data;
 37:   PetscInt       n, N, ny, Ny;

 41:   PetscLayoutGetLocalSize(is->map, &n);
 42:   PetscLayoutGetSize(is->map, &N);
 43:   PetscLayoutGetLocalSize(isy->map, &ny);
 44:   PetscLayoutGetSize(isy->map, &Ny);
 45:   if (n != ny || N != Ny) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Index sets incompatible");
 46:   PetscArraycpy(isy_general->idx,is_general->idx,n);
 47:   return(0);
 48: }

 50: static PetscErrorCode ISOnComm_General(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
 51: {
 53:   IS_General     *sub = (IS_General*)is->data;
 54:   PetscInt       n;

 57:   if (mode == PETSC_OWN_POINTER) SETERRQ(comm,PETSC_ERR_ARG_WRONG,"Cannot use PETSC_OWN_POINTER");
 58:   PetscLayoutGetLocalSize(is->map, &n);
 59:   ISCreateGeneral(comm,n,sub->idx,mode,newis);
 60:   return(0);
 61: }

 63: static PetscErrorCode ISSetBlockSize_General(IS is,PetscInt bs)
 64: {

 68:   PetscLayoutSetBlockSize(is->map, bs);
 69:   return(0);
 70: }

 72: static PetscErrorCode ISContiguousLocal_General(IS is,PetscInt gstart,PetscInt gend,PetscInt *start,PetscBool *contig)
 73: {
 74:   IS_General *sub = (IS_General*)is->data;
 75:   PetscInt   n,i,p;

 79:   *start  = 0;
 80:   *contig = PETSC_TRUE;
 81:   PetscLayoutGetLocalSize(is->map, &n);
 82:   if (!n) return(0);
 83:   p = sub->idx[0];
 84:   if (p < gstart) goto nomatch;
 85:   *start = p - gstart;
 86:   if (n > gend-p) goto nomatch;
 87:   for (i=1; i<n; i++,p++) {
 88:     if (sub->idx[i] != p+1) goto nomatch;
 89:   }
 90:   return(0);
 91: nomatch:
 92:   *start  = -1;
 93:   *contig = PETSC_FALSE;
 94:   return(0);
 95: }

 97: static PetscErrorCode ISLocate_General(IS is,PetscInt key,PetscInt *location)
 98: {
 99:   IS_General     *sub = (IS_General*)is->data;
100:   PetscInt       numIdx, i;
101:   PetscBool      sorted;

105:   PetscLayoutGetLocalSize(is->map,&numIdx);
106:   ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,&sorted);
107:   if (sorted) {PetscFindInt(key,numIdx,sub->idx,location);}
108:   else {
109:     const PetscInt *idx = sub->idx;

111:     *location = -1;
112:     for (i = 0; i < numIdx; i++) {
113:       if (idx[i] == key) {
114:         *location = i;
115:         return(0);
116:       }
117:     }
118:   }
119:   return(0);
120: }

122: static PetscErrorCode ISGetIndices_General(IS in,const PetscInt *idx[])
123: {
124:   IS_General *sub = (IS_General*)in->data;

127:   *idx = sub->idx;
128:   return(0);
129: }

131: static PetscErrorCode ISRestoreIndices_General(IS in,const PetscInt *idx[])
132: {
133:   IS_General *sub = (IS_General*)in->data;

136:   if (*idx != sub->idx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Must restore with value from ISGetIndices()");
137:   return(0);
138: }

140: static PetscErrorCode ISInvertPermutation_General(IS is,PetscInt nlocal,IS *isout)
141: {
142:   IS_General     *sub = (IS_General*)is->data;
143:   PetscInt       i,*ii,n,nstart;
144:   const PetscInt *idx = sub->idx;
145:   PetscMPIInt    size;
146:   IS             istmp,nistmp;

150:   PetscLayoutGetLocalSize(is->map, &n);
151:   MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);
152:   if (size == 1) {
153:     PetscMalloc1(n,&ii);
154:     for (i=0; i<n; i++) ii[idx[i]] = i;
155:     ISCreateGeneral(PETSC_COMM_SELF,n,ii,PETSC_OWN_POINTER,isout);
156:     ISSetPermutation(*isout);
157:   } else {
158:     /* crude, nonscalable get entire IS on each processor */
159:     ISAllGather(is,&istmp);
160:     ISSetPermutation(istmp);
161:     ISInvertPermutation(istmp,PETSC_DECIDE,&nistmp);
162:     ISDestroy(&istmp);
163:     /* get the part we need */
164:     if (nlocal == PETSC_DECIDE) nlocal = n;
165:     MPI_Scan(&nlocal,&nstart,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)is));
166: #if defined(PETSC_USE_DEBUG)
167:     {
168:       PetscInt    N;
169:       PetscMPIInt rank;
170:       MPI_Comm_rank(PetscObjectComm((PetscObject)is),&rank);
171:       PetscLayoutGetSize(is->map, &N);
172:       if (rank == size-1) {
173:         if (nstart != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Sum of nlocal lengths %d != total IS length %d",nstart,N);
174:       }
175:     }
176: #endif
177:     nstart -= nlocal;
178:     ISGetIndices(nistmp,&idx);
179:     ISCreateGeneral(PetscObjectComm((PetscObject)is),nlocal,idx+nstart,PETSC_COPY_VALUES,isout);
180:     ISRestoreIndices(nistmp,&idx);
181:     ISDestroy(&nistmp);
182:   }
183:   return(0);
184: }

186: #if defined(PETSC_HAVE_HDF5)
187: static PetscErrorCode ISView_General_HDF5(IS is, PetscViewer viewer)
188: {
189:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;
190:   hid_t           filespace;  /* file dataspace identifier */
191:   hid_t           chunkspace; /* chunk dataset property identifier */
192:   hid_t           dset_id;    /* dataset identifier */
193:   hid_t           memspace;   /* memory dataspace identifier */
194:   hid_t           inttype;    /* int type (H5T_NATIVE_INT or H5T_NATIVE_LLONG) */
195:   hid_t           file_id, group;
196:   hsize_t         dim, maxDims[3], dims[3], chunkDims[3], count[3],offset[3];
197:   PetscInt        bs, N, n, timestep, low;
198:   const PetscInt *ind;
199:   const char     *isname;
200:   PetscErrorCode  ierr;

203:   ISGetBlockSize(is,&bs);
204:   bs   = PetscMax(bs, 1); /* If N = 0, bs  = 0 as well */
205:   PetscViewerHDF5OpenGroup(viewer, &file_id, &group);
206:   PetscViewerHDF5GetTimestep(viewer, &timestep);

208:   /* Create the dataspace for the dataset.
209:    *
210:    * dims - holds the current dimensions of the dataset
211:    *
212:    * maxDims - holds the maximum dimensions of the dataset (unlimited
213:    * for the number of time steps with the current dimensions for the
214:    * other dimensions; so only additional time steps can be added).
215:    *
216:    * chunkDims - holds the size of a single time step (required to
217:    * permit extending dataset).
218:    */
219:   dim = 0;
220:   if (timestep >= 0) {
221:     dims[dim]      = timestep+1;
222:     maxDims[dim]   = H5S_UNLIMITED;
223:     chunkDims[dim] = 1;
224:     ++dim;
225:   }
226:   ISGetSize(is, &N);
227:   ISGetLocalSize(is, &n);
228:   PetscHDF5IntCast(N/bs,dims + dim);

230:   maxDims[dim]   = dims[dim];
231:   chunkDims[dim] = PetscMax(1,dims[dim]);
232:   ++dim;
233:   if (bs >= 1) {
234:     dims[dim]      = bs;
235:     maxDims[dim]   = dims[dim];
236:     chunkDims[dim] = dims[dim];
237:     ++dim;
238:   }
239:   PetscStackCallHDF5Return(filespace,H5Screate_simple,(dim, dims, maxDims));

241: #if defined(PETSC_USE_64BIT_INDICES)
242:   inttype = H5T_NATIVE_LLONG;
243: #else
244:   inttype = H5T_NATIVE_INT;
245: #endif

247:   /* Create the dataset with default properties and close filespace */
248:   PetscObjectGetName((PetscObject) is, &isname);
249:   if (!H5Lexists(group, isname, H5P_DEFAULT)) {
250:     /* Create chunk */
251:     PetscStackCallHDF5Return(chunkspace,H5Pcreate,(H5P_DATASET_CREATE));
252:     PetscStackCallHDF5(H5Pset_chunk,(chunkspace, dim, chunkDims));

254:     PetscStackCallHDF5Return(dset_id,H5Dcreate2,(group, isname, inttype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT));
255:     PetscStackCallHDF5(H5Pclose,(chunkspace));
256:   } else {
257:     PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, isname, H5P_DEFAULT));
258:     PetscStackCallHDF5(H5Dset_extent,(dset_id, dims));
259:   }
260:   PetscStackCallHDF5(H5Sclose,(filespace));

262:   /* Each process defines a dataset and writes it to the hyperslab in the file */
263:   dim = 0;
264:   if (timestep >= 0) {
265:     count[dim] = 1;
266:     ++dim;
267:   }
268:   PetscHDF5IntCast(n/bs,count + dim);
269:   ++dim;
270:   if (bs >= 1) {
271:     count[dim] = bs;
272:     ++dim;
273:   }
274:   if (n > 0) {
275:     PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL));
276:   } else {
277:     /* Can't create dataspace with zero for any dimension, so create null dataspace. */
278:     PetscStackCallHDF5Return(memspace,H5Screate,(H5S_NULL));
279:   }

281:   /* Select hyperslab in the file */
282:   PetscLayoutGetRange(is->map, &low, NULL);
283:   dim  = 0;
284:   if (timestep >= 0) {
285:     offset[dim] = timestep;
286:     ++dim;
287:   }
288:   PetscHDF5IntCast(low/bs,offset + dim);
289:   ++dim;
290:   if (bs >= 1) {
291:     offset[dim] = 0;
292:     ++dim;
293:   }
294:   if (n > 0) {
295:     PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id));
296:     PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL));
297:   } else {
298:     /* Create null filespace to match null memspace. */
299:     PetscStackCallHDF5Return(filespace,H5Screate,(H5S_NULL));
300:   }

302:   ISGetIndices(is, &ind);
303:   PetscStackCallHDF5(H5Dwrite,(dset_id, inttype, memspace, filespace, hdf5->dxpl_id, ind));
304:   PetscStackCallHDF5(H5Fflush,(file_id, H5F_SCOPE_GLOBAL));
305:   ISRestoreIndices(is, &ind);

307:   /* Close/release resources */
308:   PetscStackCallHDF5(H5Gclose,(group));
309:   PetscStackCallHDF5(H5Sclose,(filespace));
310:   PetscStackCallHDF5(H5Sclose,(memspace));
311:   PetscStackCallHDF5(H5Dclose,(dset_id));
312:   PetscInfo1(is, "Wrote IS object with name %s\n", isname);
313:   return(0);
314: }
315: #endif

317: static PetscErrorCode ISView_General_Binary(IS is,PetscViewer viewer)
318: {
320:   PetscBool      skipHeader,useMPIIO;
321:   IS_General     *isa = (IS_General*) is->data;
322:   PetscMPIInt    rank,size,mesgsize,tag = ((PetscObject)viewer)->tag, mesglen;
323:   PetscInt       n,N,len,j,tr[2];
324:   int            fdes;
325:   MPI_Status     status;
326:   PetscInt       message_count,flowcontrolcount,*values;

329:   /* ISGetLayout(is,&map); */
330:   PetscLayoutGetLocalSize(is->map, &n);
331:   PetscLayoutGetSize(is->map, &N);

333:   tr[0] = IS_FILE_CLASSID;
334:   tr[1] = N;

336:   PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
337:   if (!skipHeader) {
338:     PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT,PETSC_FALSE);
339:   }

341:   PetscViewerBinaryGetUseMPIIO(viewer,&useMPIIO);
342: #if defined(PETSC_HAVE_MPIIO)
343:   if (useMPIIO) {
344:     MPI_File       mfdes;
345:     MPI_Offset     off;
346:     PetscMPIInt    lsize;
347:     PetscInt       rstart;
348:     const PetscInt *iarray;

350:     PetscMPIIntCast(n,&lsize);
351:     PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);
352:     PetscViewerBinaryGetMPIIOOffset(viewer,&off);
353:     PetscLayoutGetRange(is->map,&rstart,NULL);
354:     off += rstart*(MPI_Offset)sizeof(PetscInt); /* off is MPI_Offset, not PetscMPIInt */
355:     MPI_File_set_view(mfdes,off,MPIU_INT,MPIU_INT,(char*)"native",MPI_INFO_NULL);
356:     ISGetIndices(is,&iarray);
357:     MPIU_File_write_all(mfdes,(void*)iarray,lsize,MPIU_INT,MPI_STATUS_IGNORE);
358:     ISRestoreIndices(is,&iarray);
359:     PetscViewerBinaryAddMPIIOOffset(viewer,N*(MPI_Offset)sizeof(PetscInt));
360:     return(0);
361:   }
362: #endif

364:   PetscViewerBinaryGetDescriptor(viewer,&fdes);
365:   MPI_Comm_rank(PetscObjectComm((PetscObject)is),&rank);
366:   MPI_Comm_size(PetscObjectComm((PetscObject)is),&size);

368:   /* determine maximum message to arrive */
369:   MPI_Reduce(&n,&len,1,MPIU_INT,MPI_MAX,0,PetscObjectComm((PetscObject)is));

371:   PetscViewerFlowControlStart(viewer,&message_count,&flowcontrolcount);
372:   if (!rank) {
373:     PetscBinaryWrite(fdes,isa->idx,n,PETSC_INT,PETSC_FALSE);

375:     PetscMalloc1(len,&values);
376:     PetscMPIIntCast(len,&mesgsize);
377:     /* receive and save messages */
378:     for (j=1; j<size; j++) {
379:       PetscViewerFlowControlStepMaster(viewer,j,&message_count,flowcontrolcount);
380:       MPI_Recv(values,mesgsize,MPIU_INT,j,tag,PetscObjectComm((PetscObject)is),&status);
381:       MPI_Get_count(&status,MPIU_INT,&mesglen);
382:       PetscBinaryWrite(fdes,values,(PetscInt)mesglen,PETSC_INT,PETSC_TRUE);
383:     }
384:     PetscViewerFlowControlEndMaster(viewer,&message_count);
385:     PetscFree(values);
386:   } else {
387:     PetscViewerFlowControlStepWorker(viewer,rank,&message_count);
388:     PetscMPIIntCast(n,&mesgsize);
389:     MPI_Send(isa->idx,mesgsize,MPIU_INT,0,tag,PetscObjectComm((PetscObject)is));
390:     PetscViewerFlowControlEndWorker(viewer,&message_count);
391:   }
392:   return(0);
393: }

395: static PetscErrorCode ISView_General(IS is,PetscViewer viewer)
396: {
397:   IS_General     *sub = (IS_General*)is->data;
399:   PetscInt       i,n,*idx = sub->idx;
400:   PetscBool      iascii,isbinary,ishdf5;

403:   PetscLayoutGetLocalSize(is->map, &n);
404:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
405:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
406:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
407:   if (iascii) {
408:     MPI_Comm          comm;
409:     PetscMPIInt       rank,size;
410:     PetscViewerFormat fmt;
411:     PetscBool         isperm;

413:     PetscObjectGetComm((PetscObject)viewer,&comm);
414:     MPI_Comm_rank(comm,&rank);
415:     MPI_Comm_size(comm,&size);

417:     PetscViewerGetFormat(viewer,&fmt);
418:     ISGetInfo(is,IS_PERMUTATION,IS_LOCAL,PETSC_FALSE,&isperm);
419:     if (isperm && fmt != PETSC_VIEWER_ASCII_MATLAB) {PetscViewerASCIIPrintf(viewer,"Index set is permutation\n");}
420:     PetscViewerASCIIPushSynchronized(viewer);
421:     if (size > 1) {
422:       if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
423:         const char* name;

425:         PetscObjectGetName((PetscObject)is,&name);
426:         PetscViewerASCIISynchronizedPrintf(viewer,"%s_%d = [...\n",name,rank);
427:         for (i=0; i<n; i++) {
428:           PetscViewerASCIISynchronizedPrintf(viewer,"%D\n",idx[i]+1);
429:         }
430:         PetscViewerASCIISynchronizedPrintf(viewer,"];\n");
431:       } else {
432:         PetscInt  st = 0;

434:         if (fmt == PETSC_VIEWER_ASCII_INDEX) st = is->map->rstart;
435:         PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in set %D\n",rank,n);
436:         for (i=0; i<n; i++) {
437:           PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %D %D\n",rank,i + st,idx[i]);
438:         }
439:       }
440:     } else {
441:       if (fmt == PETSC_VIEWER_ASCII_MATLAB) {
442:         const char* name;

444:         PetscObjectGetName((PetscObject)is,&name);
445:         PetscViewerASCIISynchronizedPrintf(viewer,"%s = [...\n",name);
446:         for (i=0; i<n; i++) {
447:           PetscViewerASCIISynchronizedPrintf(viewer,"%D\n",idx[i]+1);
448:         }
449:         PetscViewerASCIISynchronizedPrintf(viewer,"];\n");
450:       } else {
451:         PetscInt  st = 0;

453:         if (fmt == PETSC_VIEWER_ASCII_INDEX) st = is->map->rstart;
454:         PetscViewerASCIISynchronizedPrintf(viewer,"Number of indices in set %D\n",n);
455:         for (i=0; i<n; i++) {
456:           PetscViewerASCIISynchronizedPrintf(viewer,"%D %D\n",i + st,idx[i]);
457:         }
458:       }
459:     }
460:     PetscViewerFlush(viewer);
461:     PetscViewerASCIIPopSynchronized(viewer);
462:   } else if (isbinary) {
463:     ISView_General_Binary(is,viewer);
464:   } else if (ishdf5) {
465: #if defined(PETSC_HAVE_HDF5)
466:     ISView_General_HDF5(is,viewer);
467: #endif
468:   }
469:   return(0);
470: }

472: static PetscErrorCode ISSort_General(IS is)
473: {
474:   IS_General     *sub = (IS_General*)is->data;
475:   PetscInt       n;

479:   PetscLayoutGetLocalSize(is->map, &n);
480:   PetscSortInt(n,sub->idx);
481:   return(0);
482: }

484: static PetscErrorCode ISSortRemoveDups_General(IS is)
485: {
486:   IS_General     *sub = (IS_General*)is->data;
487:   PetscLayout    map;
488:   PetscInt       n;
489:   PetscBool      sorted;

493:   PetscLayoutGetLocalSize(is->map, &n);
494:   ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,&sorted);
495:   if (sorted) {
496:     PetscSortedRemoveDupsInt(&n,sub->idx);
497:   } else {
498:     PetscSortRemoveDupsInt(&n,sub->idx);
499:   }
500:   PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is), n, PETSC_DECIDE, is->map->bs, &map);
501:   PetscLayoutDestroy(&is->map);
502:   is->map = map;
503:   return(0);
504: }

506: static PetscErrorCode ISSorted_General(IS is,PetscBool  *flg)
507: {

511:   ISGetInfo(is,IS_SORTED,IS_LOCAL,PETSC_TRUE,flg);
512:   return(0);
513: }

515: PetscErrorCode  ISToGeneral_General(IS is)
516: {
518:   return(0);
519: }

521: static struct _ISOps myops = { ISGetIndices_General,
522:                                ISRestoreIndices_General,
523:                                ISInvertPermutation_General,
524:                                ISSort_General,
525:                                ISSortRemoveDups_General,
526:                                ISSorted_General,
527:                                ISDuplicate_General,
528:                                ISDestroy_General,
529:                                ISView_General,
530:                                ISLoad_Default,
531:                                ISCopy_General,
532:                                ISToGeneral_General,
533:                                ISOnComm_General,
534:                                ISSetBlockSize_General,
535:                                ISContiguousLocal_General,
536:                                ISLocate_General,
537:                                /* no specializations of {sorted,unique,perm,interval}{local,global}
538:                                 * because the default checks in ISGetInfo_XXX in index.c are exactly
539:                                 * what we would do for ISGeneral */
540:                                NULL,
541:                                NULL,
542:                                NULL,
543:                                NULL,
544:                                NULL,
545:                                NULL,
546:                                NULL,
547:                                NULL};

549: PETSC_INTERN PetscErrorCode ISSetUp_General(IS);

551: PetscErrorCode ISSetUp_General(IS is)
552: {
554:   IS_General     *sub = (IS_General*)is->data;
555:   const PetscInt *idx = sub->idx;
556:   PetscInt       n,i,min,max;

559:   PetscLayoutGetLocalSize(is->map, &n);

561:   if (n) {
562:     min = max = idx[0];
563:     for (i=1; i<n; i++) {
564:       if (idx[i] < min) min = idx[i];
565:       if (idx[i] > max) max = idx[i];
566:     }
567:     is->min = min;
568:     is->max = max;
569:   } else {
570:     is->min = PETSC_MAX_INT;
571:     is->max = PETSC_MIN_INT;
572:   }
573:   return(0);
574: }

576: /*@
577:    ISCreateGeneral - Creates a data structure for an index set
578:    containing a list of integers.

580:    Collective

582:    Input Parameters:
583: +  comm - the MPI communicator
584: .  n - the length of the index set
585: .  idx - the list of integers
586: -  mode - PETSC_COPY_VALUES, PETSC_OWN_POINTER, or PETSC_USE_POINTER; see PetscCopyMode for meaning of this flag.

588:    Output Parameter:
589: .  is - the new index set

591:    Notes:
592:    When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
593:    conceptually the same as MPI_Group operations. The IS are then
594:    distributed sets of indices and thus certain operations on them are
595:    collective.


598:    Level: beginner


601: .seealso: ISCreateStride(), ISCreateBlock(), ISAllGather(), PETSC_COPY_VALUES, PETSC_OWN_POINTER, PETSC_USE_POINTER, PetscCopyMode
602: @*/
603: PetscErrorCode  ISCreateGeneral(MPI_Comm comm,PetscInt n,const PetscInt idx[],PetscCopyMode mode,IS *is)
604: {

608:   ISCreate(comm,is);
609:   ISSetType(*is,ISGENERAL);
610:   ISGeneralSetIndices(*is,n,idx,mode);
611:   return(0);
612: }

614: /*@
615:    ISGeneralSetIndices - Sets the indices for an ISGENERAL index set

617:    Collective on IS

619:    Input Parameters:
620: +  is - the index set
621: .  n - the length of the index set
622: .  idx - the list of integers
623: -  mode - see PetscCopyMode for meaning of this flag.

625:    Level: beginner


628: .seealso: ISCreateGeneral(), ISCreateStride(), ISCreateBlock(), ISAllGather()
629: @*/
630: PetscErrorCode  ISGeneralSetIndices(IS is,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
631: {

635:   ISClearInfoCache(is,PETSC_FALSE);
636:   PetscUseMethod(is,"ISGeneralSetIndices_C",(IS,PetscInt,const PetscInt[],PetscCopyMode),(is,n,idx,mode));
637:   return(0);
638: }

640: PetscErrorCode  ISGeneralSetIndices_General(IS is,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
641: {
642:   PetscLayout    map;
644:   IS_General     *sub = (IS_General*)is->data;

647:   if (n < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"length < 0");

650:   PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)is),n,PETSC_DECIDE,is->map->bs,&map);
651:   PetscLayoutDestroy(&is->map);
652:   is->map = map;

654:   if (sub->allocated) {PetscFree(sub->idx);}
655:   if (mode == PETSC_COPY_VALUES) {
656:     PetscMalloc1(n,&sub->idx);
657:     PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
658:     PetscArraycpy(sub->idx,idx,n);
659:     sub->allocated = PETSC_TRUE;
660:   } else if (mode == PETSC_OWN_POINTER) {
661:     sub->idx = (PetscInt*)idx;
662:     PetscLogObjectMemory((PetscObject)is,n*sizeof(PetscInt));
663:     sub->allocated = PETSC_TRUE;
664:   } else {
665:     sub->idx = (PetscInt*)idx;
666:     sub->allocated = PETSC_FALSE;
667:   }

669:   ISSetUp_General(is);
670:   ISViewFromOptions(is,NULL,"-is_view");
671:   return(0);
672: }

674: static PetscErrorCode ISGeneralFilter_General(IS is, PetscInt start, PetscInt end)
675: {
676:   IS_General     *sub = (IS_General*)is->data;
677:   PetscInt       *idx = sub->idx,*idxnew;
678:   PetscInt       i,n = is->map->n,nnew = 0,o;

682:   for (i=0; i<n; ++i)
683:     if (idx[i] >= start && idx[i] < end)
684:       nnew++;
685:   PetscMalloc1(nnew, &idxnew);
686:   for (o=0, i=0; i<n; i++) {
687:     if (idx[i] >= start && idx[i] < end)
688:       idxnew[o++] = idx[i];
689:   }
690:   ISGeneralSetIndices_General(is,nnew,idxnew,PETSC_OWN_POINTER);
691:   return(0);
692: }

694: /*@
695:    ISGeneralFilter - Remove all points outside of [start, end)

697:    Collective on IS

699:    Input Parameters:
700: +  is - the index set
701: .  start - the lowest index kept
702: -  end - one more than the highest index kept

704:    Level: beginner

706: .seealso: ISCreateGeneral(), ISGeneralSetIndices()
707: @*/
708: PetscErrorCode ISGeneralFilter(IS is, PetscInt start, PetscInt end)
709: {

714:   ISClearInfoCache(is,PETSC_FALSE);
715:   PetscUseMethod(is,"ISGeneralFilter_C",(IS,PetscInt,PetscInt),(is,start,end));
716:   return(0);
717: }

719: PETSC_EXTERN PetscErrorCode ISCreate_General(IS is)
720: {
722:   IS_General     *sub;

725:   PetscNewLog(is,&sub);
726:   is->data = (void *) sub;
727:   PetscMemcpy(is->ops,&myops,sizeof(myops));
728:   PetscObjectComposeFunction((PetscObject)is,"ISGeneralSetIndices_C",ISGeneralSetIndices_General);
729:   PetscObjectComposeFunction((PetscObject)is,"ISGeneralFilter_C",ISGeneralFilter_General);
730:   return(0);
731: }