Actual source code: block.c
petsc-3.3-p7 2013-05-11
2: /*
3: Provides the functions for index sets (IS) defined by a list of integers.
4: These are for blocks of data, each block is indicated with a single integer.
5: */
6: #include <petsc-private/isimpl.h> /*I "petscis.h" I*/
7: #include <petscvec.h>
9: typedef struct {
10: PetscInt N,n; /* number of blocks */
11: PetscBool sorted; /* are the blocks sorted? */
12: PetscInt *idx;
13: } IS_Block;
17: PetscErrorCode ISDestroy_Block(IS is)
18: {
19: IS_Block *is_block = (IS_Block*)is->data;
23: PetscFree(is_block->idx);
24: PetscObjectComposeFunctionDynamic((PetscObject)is,"ISBlockSetIndices_C","",0);
25: PetscObjectComposeFunctionDynamic((PetscObject)is,"ISBlockGetIndices_C","",0);
26: PetscObjectComposeFunctionDynamic((PetscObject)is,"ISBlockRestoreIndices_C","",0);
27: PetscObjectComposeFunctionDynamic((PetscObject)is,"ISBlockGetSize_C","",0);
28: PetscObjectComposeFunctionDynamic((PetscObject)is,"ISBlockGetLocalSize_C","",0);
29: PetscFree(is->data);
30: return(0);
31: }
35: PetscErrorCode ISGetIndices_Block(IS in,const PetscInt *idx[])
36: {
37: IS_Block *sub = (IS_Block*)in->data;
39: PetscInt i,j,k,bs = in->bs,n = sub->n,*ii,*jj;
42: if (bs == 1) {
43: *idx = sub->idx;
44: } else {
45: PetscMalloc(bs*n*sizeof(PetscInt),&jj);
46: *idx = jj;
47: k = 0;
48: ii = sub->idx;
49: for (i=0; i<n; i++) {
50: for (j=0; j<bs; j++) {
51: jj[k++] = bs*ii[i] + j;
52: }
53: }
54: }
55: return(0);
56: }
60: PetscErrorCode ISRestoreIndices_Block(IS in,const PetscInt *idx[])
61: {
62: IS_Block *sub = (IS_Block*)in->data;
66: if (in->bs != 1) {
67: PetscFree(*(void**)idx);
68: } else {
69: if (*idx != sub->idx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Must restore with value from ISGetIndices()");
70: }
71: return(0);
72: }
76: PetscErrorCode ISGetSize_Block(IS is,PetscInt *size)
77: {
78: IS_Block *sub = (IS_Block *)is->data;
81: *size = is->bs*sub->N;
82: return(0);
83: }
87: PetscErrorCode ISGetLocalSize_Block(IS is,PetscInt *size)
88: {
89: IS_Block *sub = (IS_Block *)is->data;
92: *size = is->bs*sub->n;
93: return(0);
94: }
98: PetscErrorCode ISInvertPermutation_Block(IS is,PetscInt nlocal,IS *isout)
99: {
100: IS_Block *sub = (IS_Block *)is->data;
101: PetscInt i,*ii,n = sub->n,*idx = sub->idx;
102: PetscMPIInt size;
106: MPI_Comm_size(((PetscObject)is)->comm,&size);
107: if (size == 1) {
108: PetscMalloc(n*sizeof(PetscInt),&ii);
109: for (i=0; i<n; i++) {
110: ii[idx[i]] = i;
111: }
112: ISCreateBlock(PETSC_COMM_SELF,is->bs,n,ii,PETSC_OWN_POINTER,isout);
113: ISSetPermutation(*isout);
114: } else {
115: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No inversion written yet for block IS");
116: }
117: return(0);
118: }
122: PetscErrorCode ISView_Block(IS is, PetscViewer viewer)
123: {
124: IS_Block *sub = (IS_Block *)is->data;
126: PetscInt i,n = sub->n,*idx = sub->idx;
127: PetscBool iascii;
130: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
131: if (iascii) {
132: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
133: if (is->isperm) {
134: PetscViewerASCIISynchronizedPrintf(viewer,"Block Index set is permutation\n");
135: }
136: PetscViewerASCIISynchronizedPrintf(viewer,"Block size %D\n",is->bs);
137: PetscViewerASCIISynchronizedPrintf(viewer,"Number of block indices in set %D\n",n);
138: PetscViewerASCIISynchronizedPrintf(viewer,"The first indices of each block are\n");
139: for (i=0; i<n; i++) {
140: PetscViewerASCIISynchronizedPrintf(viewer,"Block %D Index %D\n",i,idx[i]);
141: }
142: PetscViewerFlush(viewer);
143: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
144: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for this object",((PetscObject)viewer)->type_name);
145: return(0);
146: }
150: PetscErrorCode ISSort_Block(IS is)
151: {
152: IS_Block *sub = (IS_Block *)is->data;
156: if (sub->sorted) return(0);
157: PetscSortInt(sub->n,sub->idx);
158: sub->sorted = PETSC_TRUE;
159: return(0);
160: }
164: PetscErrorCode ISSorted_Block(IS is,PetscBool *flg)
165: {
166: IS_Block *sub = (IS_Block *)is->data;
169: *flg = sub->sorted;
170: return(0);
171: }
175: PetscErrorCode ISDuplicate_Block(IS is,IS *newIS)
176: {
178: IS_Block *sub = (IS_Block *)is->data;
181: ISCreateBlock(((PetscObject)is)->comm,is->bs,sub->n,sub->idx,PETSC_COPY_VALUES,newIS);
182: return(0);
183: }
187: PetscErrorCode ISIdentity_Block(IS is,PetscBool *ident)
188: {
189: IS_Block *is_block = (IS_Block*)is->data;
190: PetscInt i,n = is_block->n,*idx = is_block->idx,bs = is->bs;
193: is->isidentity = PETSC_TRUE;
194: *ident = PETSC_TRUE;
195: for (i=0; i<n; i++) {
196: if (idx[i] != bs*i) {
197: is->isidentity = PETSC_FALSE;
198: *ident = PETSC_FALSE;
199: return(0);
200: }
201: }
202: return(0);
203: }
207: static PetscErrorCode ISCopy_Block(IS is,IS isy)
208: {
209: IS_Block *is_block = (IS_Block*)is->data,*isy_block = (IS_Block*)isy->data;
213: if (is_block->n != isy_block->n || is_block->N != isy_block->N || is->bs != isy->bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Index sets incompatible");
214: isy_block->sorted = is_block->sorted;
215: PetscMemcpy(isy_block->idx,is_block->idx,is_block->n*sizeof(PetscInt));
216: return(0);
217: }
221: static PetscErrorCode ISOnComm_Block(IS is,MPI_Comm comm,PetscCopyMode mode,IS *newis)
222: {
224: IS_Block *sub = (IS_Block*)is->data;
227: if (mode == PETSC_OWN_POINTER) SETERRQ(comm,PETSC_ERR_ARG_WRONG,"Cannot use PETSC_OWN_POINTER");
228: ISCreateBlock(comm,is->bs,sub->n,sub->idx,mode,newis);
229: return(0);
230: }
234: static PetscErrorCode ISSetBlockSize_Block(IS is,PetscInt bs)
235: {
238: if (is->bs != bs) SETERRQ2(((PetscObject)is)->comm,PETSC_ERR_ARG_SIZ,"Cannot change block size for ISBLOCK from %D to %D",is->bs,bs);
239: return(0);
240: }
243: static struct _ISOps myops = { ISGetSize_Block,
244: ISGetLocalSize_Block,
245: ISGetIndices_Block,
246: ISRestoreIndices_Block,
247: ISInvertPermutation_Block,
248: ISSort_Block,
249: ISSorted_Block,
250: ISDuplicate_Block,
251: ISDestroy_Block,
252: ISView_Block,
253: ISIdentity_Block,
254: ISCopy_Block,
255: 0,
256: ISOnComm_Block,
257: ISSetBlockSize_Block
258: };
262: /*@
263: ISBlockSetIndices - The indices are relative to entries, not blocks.
265: Collective on IS
267: Input Parameters:
268: + is - the index set
269: . bs - number of elements in each block, one for each block and count of block not indices
270: . n - the length of the index set (the number of blocks)
271: . idx - the list of integers, these are by block, not by location
272: + mode - see PetscCopyMode, only PETSC_COPY_VALUES and PETSC_OWN_POINTER are supported
275: Notes:
276: When the communicator is not MPI_COMM_SELF, the operations on the
277: index sets, IS, are NOT conceptually the same as MPI_Group operations.
278: The index sets are then distributed sets of indices and thus certain operations
279: on them are collective.
281: Example:
282: If you wish to index the values {0,1,4,5}, then use
283: a block size of 2 and idx of {0,2}.
285: Level: beginner
287: Concepts: IS^block
288: Concepts: index sets^block
289: Concepts: block^index set
291: .seealso: ISCreateStride(), ISCreateGeneral(), ISAllGather()
292: @*/
293: PetscErrorCode ISBlockSetIndices(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
294: {
297: PetscUseMethod(is,"ISBlockSetIndices_C",(IS,PetscInt,PetscInt,const PetscInt[],PetscCopyMode),(is,bs,n,idx,mode));
298: return(0);
299: }
301: EXTERN_C_BEGIN
304: PetscErrorCode ISBlockSetIndices_Block(IS is,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode)
305: {
307: PetscInt i,min,max;
308: IS_Block *sub = (IS_Block*)is->data;
309: PetscBool sorted = PETSC_TRUE;
312: PetscFree(sub->idx);
313: sub->n = n;
314: MPI_Allreduce(&n,&sub->N,1,MPIU_INT,MPI_SUM,((PetscObject)is)->comm);
315: for (i=1; i<n; i++) {
316: if (idx[i] < idx[i-1]) {sorted = PETSC_FALSE; break;}
317: }
318: if (n) {min = max = idx[0];} else {min = max = 0;}
319: for (i=1; i<n; i++) {
320: if (idx[i] < min) min = idx[i];
321: if (idx[i] > max) max = idx[i];
322: }
323: if (mode == PETSC_COPY_VALUES) {
324: PetscMalloc(n*sizeof(PetscInt),&sub->idx);
325: PetscLogObjectMemory(is,n*sizeof(PetscInt));
326: PetscMemcpy(sub->idx,idx,n*sizeof(PetscInt));
327: } else if (mode == PETSC_OWN_POINTER) {
328: sub->idx = (PetscInt*) idx;
329: } else SETERRQ(((PetscObject)is)->comm,PETSC_ERR_SUP,"Only supports PETSC_COPY_VALUES and PETSC_OWN_POINTER");
330: sub->sorted = sorted;
331: is->bs = bs;
332: is->min = bs*min;
333: is->max = bs*max+bs-1;
334: is->data = (void*)sub;
335: PetscMemcpy(is->ops,&myops,sizeof(myops));
336: is->isperm = PETSC_FALSE;
337: return(0);
338: }
339: EXTERN_C_END
343: /*@
344: ISCreateBlock - Creates a data structure for an index set containing
345: a list of integers. The indices are relative to entries, not blocks.
347: Collective on MPI_Comm
349: Input Parameters:
350: + comm - the MPI communicator
351: . bs - number of elements in each block
352: . n - the length of the index set (the number of blocks)
353: . idx - the list of integers, one for each block and count of block not indices
354: - mode - see PetscCopyMode, only PETSC_COPY_VALUES and PETSC_OWN_POINTER are supported in this routine
356: Output Parameter:
357: . is - the new index set
359: Notes:
360: When the communicator is not MPI_COMM_SELF, the operations on the
361: index sets, IS, are NOT conceptually the same as MPI_Group operations.
362: The index sets are then distributed sets of indices and thus certain operations
363: on them are collective.
365: Example:
366: If you wish to index the values {0,1,6,7}, then use
367: a block size of 2 and idx of {0,3}.
369: Level: beginner
371: Concepts: IS^block
372: Concepts: index sets^block
373: Concepts: block^index set
375: .seealso: ISCreateStride(), ISCreateGeneral(), ISAllGather()
376: @*/
377: PetscErrorCode ISCreateBlock(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscInt idx[],PetscCopyMode mode,IS *is)
378: {
383: if (n < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"length < 0");
386: ISCreate(comm,is);
387: ISSetType(*is,ISBLOCK);
388: ISBlockSetIndices(*is,bs,n,idx,mode);
389: return(0);
390: }
392: EXTERN_C_BEGIN
395: PetscErrorCode ISBlockGetIndices_Block(IS is,const PetscInt *idx[])
396: {
397: IS_Block *sub = (IS_Block*)is->data;
400: *idx = sub->idx;
401: return(0);
402: }
403: EXTERN_C_END
405: EXTERN_C_BEGIN
408: PetscErrorCode ISBlockRestoreIndices_Block(IS is,const PetscInt *idx[])
409: {
411: return(0);
412: }
413: EXTERN_C_END
417: /*@C
418: ISBlockGetIndices - Gets the indices associated with each block.
420: Not Collective
422: Input Parameter:
423: . is - the index set
425: Output Parameter:
426: . idx - the integer indices, one for each block and count of block not indices
428: Level: intermediate
430: Concepts: IS^block
431: Concepts: index sets^getting indices
432: Concepts: index sets^block
434: .seealso: ISGetIndices(), ISBlockRestoreIndices()
435: @*/
436: PetscErrorCode ISBlockGetIndices(IS is,const PetscInt *idx[])
437: {
440: PetscUseMethod(is,"ISBlockGetIndices_C",(IS,const PetscInt*[]),(is,idx));
441: return(0);
442: }
446: /*@C
447: ISBlockRestoreIndices - Restores the indices associated with each block.
449: Not Collective
451: Input Parameter:
452: . is - the index set
454: Output Parameter:
455: . idx - the integer indices
457: Level: intermediate
459: Concepts: IS^block
460: Concepts: index sets^getting indices
461: Concepts: index sets^block
463: .seealso: ISRestoreIndices(), ISBlockGetIndices()
464: @*/
465: PetscErrorCode ISBlockRestoreIndices(IS is,const PetscInt *idx[])
466: {
469: PetscUseMethod(is,"ISBlockRestoreIndices_C",(IS,const PetscInt*[]),(is,idx));
470: return(0);
471: }
475: /*@
476: ISBlockGetLocalSize - Returns the local number of blocks in the index set.
478: Not Collective
480: Input Parameter:
481: . is - the index set
483: Output Parameter:
484: . size - the local number of blocks
486: Level: intermediate
488: Concepts: IS^block sizes
489: Concepts: index sets^block sizes
491: .seealso: ISGetBlockSize(), ISBlockGetSize(), ISGetSize(), ISCreateBlock()
492: @*/
493: PetscErrorCode ISBlockGetLocalSize(IS is,PetscInt *size)
494: {
497: PetscUseMethod(is,"ISBlockGetLocalSize_C",(IS,PetscInt*),(is,size));
498: return(0);
499: }
501: EXTERN_C_BEGIN
504: PetscErrorCode ISBlockGetLocalSize_Block(IS is,PetscInt *size)
505: {
506: IS_Block *sub = (IS_Block *)is->data;
509: *size = sub->n;
510: return(0);
511: }
512: EXTERN_C_END
516: /*@
517: ISBlockGetSize - Returns the global number of blocks in the index set.
519: Not Collective
521: Input Parameter:
522: . is - the index set
524: Output Parameter:
525: . size - the global number of blocks
527: Level: intermediate
529: Concepts: IS^block sizes
530: Concepts: index sets^block sizes
532: .seealso: ISGetBlockSize(), ISBlockGetLocalSize(), ISGetSize(), ISCreateBlock()
533: @*/
534: PetscErrorCode ISBlockGetSize(IS is,PetscInt *size)
535: {
538: PetscUseMethod(is,"ISBlockGetSize_C",(IS,PetscInt*),(is,size));
539: return(0);
540: }
542: EXTERN_C_BEGIN
545: PetscErrorCode ISBlockGetSize_Block(IS is,PetscInt *size)
546: {
547: IS_Block *sub = (IS_Block *)is->data;
550: *size = sub->N;
551: return(0);
552: }
553: EXTERN_C_END
557: PetscErrorCode ISToGeneral_Block(IS inis)
558: {
560: const PetscInt *idx;
561: PetscInt n;
564: ISGetLocalSize(inis,&n);
565: ISGetIndices(inis,&idx);
566: ISSetType(inis,ISGENERAL);
567: ISGeneralSetIndices(inis,n,idx,PETSC_OWN_POINTER);
568: return(0);
569: }
571: EXTERN_C_BEGIN
574: PetscErrorCode ISCreate_Block(IS is)
575: {
577: IS_Block *sub;
580: PetscNewLog(is,IS_Block,&sub);
581: is->data = sub;
582: PetscObjectComposeFunctionDynamic((PetscObject)is,"ISBlockSetIndices_C","ISBlockSetIndices_Block",
583: ISBlockSetIndices_Block);
584: PetscObjectComposeFunctionDynamic((PetscObject)is,"ISBlockGetIndices_C","ISBlockGetIndices_Block",
585: ISBlockGetIndices_Block);
586: PetscObjectComposeFunctionDynamic((PetscObject)is,"ISBlockRestoreIndices_C","ISBlockRestoreIndices_Block",
587: ISBlockRestoreIndices_Block);
588: PetscObjectComposeFunctionDynamic((PetscObject)is,"ISBlockGetSize_C","ISBlockGetSize_Block",
589: ISBlockGetSize_Block);
590: PetscObjectComposeFunctionDynamic((PetscObject)is,"ISBlockGetLocalSize_C","ISBlockGetLocalSize_Block",
591: ISBlockGetLocalSize_Block);
592: return(0);
593: }
594: EXTERN_C_END