Actual source code: general.c

  1: /*$Id: general.c,v 1.103 2001/03/23 23:21:08 balay Exp $*/
  2: /*
  3:      Provides the functions for index sets (IS) defined by a list of integers.
  4: */
 5:  #include src/vec/is/impls/general/general.h

  7: EXTERN int VecInitializePackage(char *);

  9: #undef __FUNCT__  
 11: int ISDuplicate_General(IS is,IS *newIS)
 12: {
 13:   int        ierr;
 14:   IS_General *sub = (IS_General *)is->data;

 17:   ISCreateGeneral(is->comm,sub->n,sub->idx,newIS);
 18:   return(0);
 19: }

 21: #undef __FUNCT__  
 23: int ISDestroy_General(IS is)
 24: {
 25:   IS_General *is_general = (IS_General*)is->data;
 26:   int        ierr;

 29:   PetscFree(is_general->idx);
 30:   PetscFree(is_general);
 31:   PetscLogObjectDestroy(is);
 32:   PetscHeaderDestroy(is);
 33:   return(0);
 34: }

 36: #undef __FUNCT__  
 38: int ISIdentity_General(IS is,PetscTruth *ident)
 39: {
 40:   IS_General *is_general = (IS_General*)is->data;
 41:   int        i,n = is_general->n,*idx = is_general->idx;

 44:   is->isidentity = PETSC_TRUE;
 45:   *ident         = PETSC_TRUE;
 46:   for (i=0; i<n; i++) {
 47:     if (idx[i] != i) {
 48:       is->isidentity = PETSC_FALSE;
 49:       *ident         = PETSC_FALSE;
 50:       break;
 51:     }
 52:   }
 53:   return(0);
 54: }

 56: #undef __FUNCT__  
 58: int ISGetIndices_General(IS in,int **idx)
 59: {
 60:   IS_General *sub = (IS_General*)in->data;

 63:   *idx = sub->idx;
 64:   return(0);
 65: }

 67: #undef __FUNCT__  
 69: int ISRestoreIndices_General(IS in,int **idx)
 70: {
 71:   IS_General *sub = (IS_General*)in->data;

 74:   if (*idx != sub->idx) {
 75:     SETERRQ(PETSC_ERR_ARG_WRONG,"Must restore with value from ISGetIndices()");
 76:   }
 77:   return(0);
 78: }

 80: #undef __FUNCT__  
 82: int ISGetSize_General(IS is,int *size)
 83: {
 84:   IS_General *sub = (IS_General *)is->data;

 87:   *size = sub->N;
 88:   return(0);
 89: }

 91: #undef __FUNCT__  
 93: int ISGetLocalSize_General(IS is,int *size)
 94: {
 95:   IS_General *sub = (IS_General *)is->data;

 98:   *size = sub->n;
 99:   return(0);
100: }

102: #undef __FUNCT__  
104: int ISInvertPermutation_General(IS is,int nlocal,IS *isout)
105: {
106:   IS_General *sub = (IS_General *)is->data;
107:   int        i,ierr,*ii,n = sub->n,*idx = sub->idx,size,nstart;
108:   IS         istmp,nistmp;

111:   MPI_Comm_size(is->comm,&size);
112:   if (size == 1) {
113:     PetscMalloc(n*sizeof(int),&ii);
114:     for (i=0; i<n; i++) {
115:       ii[idx[i]] = i;
116:     }
117:     ISCreateGeneral(PETSC_COMM_SELF,n,ii,isout);
118:     ISSetPermutation(*isout);
119:     PetscFree(ii);
120:   } else {
121:     /* crude, nonscalable get entire IS on each processor */
122:     ISAllGather(is,&istmp);
123:     ISSetPermutation(istmp);
124:     ISInvertPermutation(istmp,PETSC_DECIDE,&nistmp);
125:     ISDestroy(istmp);
126:     /* get the part we need */
127:     ierr    = MPI_Scan(&nlocal,&nstart,1,MPI_INT,MPI_SUM,is->comm);
128:     nstart -= nlocal;
129:     ierr    = ISGetIndices(nistmp,&idx);
130:     ierr    = ISCreateGeneral(is->comm,nlocal,idx+nstart,isout);
131:     ierr    = ISRestoreIndices(nistmp,&idx);
132:     ierr    = ISDestroy(nistmp);
133:   }
134:   return(0);
135: }

137: #undef __FUNCT__  
139: int ISView_General(IS is,PetscViewer viewer)
140: {
141:   IS_General  *sub = (IS_General *)is->data;
142:   int         i,n = sub->n,*idx = sub->idx,ierr;
143:   PetscTruth  isascii;

146:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
147:   if (isascii) {
148:     MPI_Comm comm;
149:     int      rank,size;

151:     PetscObjectGetComm((PetscObject)viewer,&comm);
152:     MPI_Comm_rank(comm,&rank);
153:     MPI_Comm_size(comm,&size);

155:     if (size > 1) {
156:       if (is->isperm) {
157:         PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Index set is permutationn",rank);
158:       }
159:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Number of indices in set %dn",rank,n);
160:       for (i=0; i<n; i++) {
161:         PetscViewerASCIISynchronizedPrintf(viewer,"[%d] %d %dn",rank,i,idx[i]);
162:       }
163:     } else {
164:       if (is->isperm) {
165:         PetscViewerASCIISynchronizedPrintf(viewer,"Index set is permutationn");
166:       }
167:       PetscViewerASCIISynchronizedPrintf(viewer,"Number of indices in set %dn",n);
168:       for (i=0; i<n; i++) {
169:         PetscViewerASCIISynchronizedPrintf(viewer,"%d %dn",i,idx[i]);
170:       }
171:     }
172:     PetscViewerFlush(viewer);
173:   } else {
174:     SETERRQ1(1,"Viewer type %s not supported for this object",((PetscObject)viewer)->type_name);
175:   }
176:   return(0);
177: }

179: #undef __FUNCT__  
181: int ISSort_General(IS is)
182: {
183:   IS_General *sub = (IS_General *)is->data;
184:   int        ierr;

187:   if (sub->sorted) return(0);
188:   PetscSortInt(sub->n,sub->idx);
189:   sub->sorted = PETSC_TRUE;
190:   return(0);
191: }

193: #undef __FUNCT__  
195: int ISSorted_General(IS is,PetscTruth *flg)
196: {
197:   IS_General *sub = (IS_General *)is->data;

200:   *flg = sub->sorted;
201:   return(0);
202: }

204: static struct _ISOps myops = { ISGetSize_General,
205:                                ISGetLocalSize_General,
206:                                ISGetIndices_General,
207:                                ISRestoreIndices_General,
208:                                ISInvertPermutation_General,
209:                                ISSort_General,
210:                                ISSorted_General,
211:                                ISDuplicate_General,
212:                                ISDestroy_General,
213:                                ISView_General,
214:                                ISIdentity_General };

216: #undef __FUNCT__  
218: /*@C
219:    ISCreateGeneral - Creates a data structure for an index set 
220:    containing a list of integers.

222:    Collective on MPI_Comm

224:    Input Parameters:
225: +  comm - the MPI communicator
226: .  n - the length of the index set
227: -  idx - the list of integers

229:    Output Parameter:
230: .  is - the new index set

232:    Notes:
233:    When the communicator is not MPI_COMM_SELF, the operations on IS are NOT
234:    conceptually the same as MPI_Group operations. The IS are then
235:    distributed sets of indices and thus certain operations on them are collective.

237:    Level: beginner

239:   Concepts: index sets^creating
240:   Concepts: IS^creating

242: .seealso: ISCreateStride(), ISCreateBlock(), ISAllGather()
243: @*/
244: int ISCreateGeneral(MPI_Comm comm,int n,const int idx[],IS *is)
245: {
246:   int        i,min,max,ierr;
247:   PetscTruth sorted = PETSC_TRUE;
248:   IS         Nindex;
249:   IS_General *sub;
250:   PetscTruth flg;

254:   if (n < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
256:   *is = PETSC_NULL;
257: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
258:   VecInitializePackage(PETSC_NULL);
259: #endif

261:   PetscHeaderCreate(Nindex,_p_IS,struct _ISOps,IS_COOKIE,IS_GENERAL,"IS",comm,ISDestroy,ISView);
262:   PetscLogObjectCreate(Nindex);
263:   ierr           = PetscNew(IS_General,&sub);
264:   PetscLogObjectMemory(Nindex,sizeof(IS_General)+n*sizeof(int)+sizeof(struct _p_IS));
265:   ierr           = PetscMalloc((n+1)*sizeof(int),&sub->idx);
266:   sub->n         = n;

268:   MPI_Allreduce(&n,&sub->N,1,MPI_INT,MPI_SUM,comm);
269:   for (i=1; i<n; i++) {
270:     if (idx[i] < idx[i-1]) {sorted = PETSC_FALSE; break;}
271:   }
272:   if (n) {min = max = idx[0];} else {min = max = 0;}
273:   for (i=1; i<n; i++) {
274:     if (idx[i] < min) min = idx[i];
275:     if (idx[i] > max) max = idx[i];
276:   }
277:   PetscMemcpy(sub->idx,idx,n*sizeof(int));
278:   sub->sorted     = sorted;
279:   Nindex->min     = min;
280:   Nindex->max     = max;
281:   Nindex->data    = (void*)sub;
282:   PetscMemcpy(Nindex->ops,&myops,sizeof(myops));
283:   Nindex->isperm     = PETSC_FALSE;
284:   Nindex->isidentity = PETSC_FALSE;
285:   PetscOptionsHasName(PETSC_NULL,"-is_view",&flg);
286:   if (flg) {
287:     ISView(Nindex,PETSC_VIEWER_STDOUT_(Nindex->comm));
288:   }
289:   *is = Nindex;
290:   return(0);
291: }