Actual source code: sorti.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  2: /*
  3:    This file contains routines for sorting integers. Values are sorted in place.
  4:  */
  5: #include <petsc/private/petscimpl.h>                /*I  "petscsys.h"  I*/

  7: #define SWAP(a,b,t) {t=a;a=b;b=t;}

  9: #define MEDIAN3(v,a,b,c)                        \
 10:   (v[a]<v[b]                                    \
 11:    ? (v[b]<v[c]                                 \
 12:       ? b                                       \
 13:       : (v[a]<v[c] ? c : a))                    \
 14:    : (v[c]<v[b]                                 \
 15:       ? b                                       \
 16:       : (v[a]<v[c] ? a : c)))

 18: #define MEDIAN(v,right) MEDIAN3(v,right/4,right/2,right/4*3)

 20: /* -----------------------------------------------------------------------*/

 24: /*
 25:    A simple version of quicksort; taken from Kernighan and Ritchie, page 87.
 26:    Assumes 0 origin for v, number of elements = right+1 (right is index of
 27:    right-most member).
 28: */
 29: static void PetscSortInt_Private(PetscInt *v,PetscInt right)
 30: {
 31:   PetscInt i,j,pivot,tmp;

 33:   if (right <= 1) {
 34:     if (right == 1) {
 35:       if (v[0] > v[1]) SWAP(v[0],v[1],tmp);
 36:     }
 37:     return;
 38:   }
 39:   i = MEDIAN(v,right);          /* Choose a pivot */
 40:   SWAP(v[0],v[i],tmp);          /* Move it out of the way */
 41:   pivot = v[0];
 42:   for (i=0,j=right+1;; ) {
 43:     while (++i < j && v[i] <= pivot) ; /* Scan from the left */
 44:     while (v[--j] > pivot) ;           /* Scan from the right */
 45:     if (i >= j) break;
 46:     SWAP(v[i],v[j],tmp);
 47:   }
 48:   SWAP(v[0],v[j],tmp);          /* Put pivot back in place. */
 49:   PetscSortInt_Private(v,j-1);
 50:   PetscSortInt_Private(v+j+1,right-(j+1));
 51: }

 55: /*@
 56:    PetscSortInt - Sorts an array of integers in place in increasing order.

 58:    Not Collective

 60:    Input Parameters:
 61: +  n  - number of values
 62: -  i  - array of integers

 64:    Level: intermediate

 66:    Concepts: sorting^ints

 68: .seealso: PetscSortReal(), PetscSortIntWithPermutation()
 69: @*/
 70: PetscErrorCode  PetscSortInt(PetscInt n,PetscInt i[])
 71: {
 72:   PetscInt j,k,tmp,ik;

 75:   if (n<8) {
 76:     for (k=0; k<n; k++) {
 77:       ik = i[k];
 78:       for (j=k+1; j<n; j++) {
 79:         if (ik > i[j]) {
 80:           SWAP(i[k],i[j],tmp);
 81:           ik = i[k];
 82:         }
 83:       }
 84:     }
 85:   } else PetscSortInt_Private(i,n-1);
 86:   return(0);
 87: }

 91: /*@
 92:    PetscSortRemoveDupsInt - Sorts an array of integers in place in increasing order removes all duplicate entries

 94:    Not Collective

 96:    Input Parameters:
 97: +  n  - number of values
 98: -  ii  - array of integers

100:    Output Parameter:
101: .  n - number of non-redundant values

103:    Level: intermediate

105:    Concepts: sorting^ints

107: .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt()
108: @*/
109: PetscErrorCode  PetscSortRemoveDupsInt(PetscInt *n,PetscInt ii[])
110: {
112:   PetscInt       i,s = 0,N = *n, b = 0;

115:   PetscSortInt(N,ii);
116:   for (i=0; i<N-1; i++) {
117:     if (ii[b+s+1] != ii[b]) {
118:       ii[b+1] = ii[b+s+1]; b++;
119:     } else s++;
120:   }
121:   *n = N - s;
122:   return(0);
123: }

127: /*@
128:   PetscFindInt - Finds integer in a sorted array of integers

130:    Not Collective

132:    Input Parameters:
133: +  key - the integer to locate
134: .  n   - number of values in the array
135: -  ii  - array of integers

137:    Output Parameter:
138: .  loc - the location if found, otherwise -(slot+1) where slot is the place the value would go

140:    Level: intermediate

142:    Concepts: sorting^ints

144: .seealso: PetscSortInt(), PetscSortIntWithArray(), PetscSortRemoveDupsInt()
145: @*/
146: PetscErrorCode PetscFindInt(PetscInt key, PetscInt n, const PetscInt ii[], PetscInt *loc)
147: {
148:   PetscInt lo = 0,hi = n;

152:   if (!n) {*loc = -1; return(0);}
154:   while (hi - lo > 1) {
155:     PetscInt mid = lo + (hi - lo)/2;
156:     if (key < ii[mid]) hi = mid;
157:     else               lo = mid;
158:   }
159:   *loc = key == ii[lo] ? lo : -(lo + (key > ii[lo]) + 1);
160:   return(0);
161: }


164: /* -----------------------------------------------------------------------*/
165: #define SWAP2(a,b,c,d,t) {t=a;a=b;b=t;t=c;c=d;d=t;}

169: /*
170:    A simple version of quicksort; taken from Kernighan and Ritchie, page 87.
171:    Assumes 0 origin for v, number of elements = right+1 (right is index of
172:    right-most member).
173: */
174: static PetscErrorCode PetscSortIntWithArray_Private(PetscInt *v,PetscInt *V,PetscInt right)
175: {
177:   PetscInt       i,vl,last,tmp;

180:   if (right <= 1) {
181:     if (right == 1) {
182:       if (v[0] > v[1]) SWAP2(v[0],v[1],V[0],V[1],tmp);
183:     }
184:     return(0);
185:   }
186:   SWAP2(v[0],v[right/2],V[0],V[right/2],tmp);
187:   vl   = v[0];
188:   last = 0;
189:   for (i=1; i<=right; i++) {
190:     if (v[i] < vl) {last++; SWAP2(v[last],v[i],V[last],V[i],tmp);}
191:   }
192:   SWAP2(v[0],v[last],V[0],V[last],tmp);
193:   PetscSortIntWithArray_Private(v,V,last-1);
194:   PetscSortIntWithArray_Private(v+last+1,V+last+1,right-(last+1));
195:   return(0);
196: }

200: /*@
201:    PetscSortIntWithArray - Sorts an array of integers in place in increasing order;
202:        changes a second array to match the sorted first array.

204:    Not Collective

206:    Input Parameters:
207: +  n  - number of values
208: .  i  - array of integers
209: -  I - second array of integers

211:    Level: intermediate

213:    Concepts: sorting^ints with array

215: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt()
216: @*/
217: PetscErrorCode  PetscSortIntWithArray(PetscInt n,PetscInt i[],PetscInt Ii[])
218: {
220:   PetscInt       j,k,tmp,ik;

223:   if (n<8) {
224:     for (k=0; k<n; k++) {
225:       ik = i[k];
226:       for (j=k+1; j<n; j++) {
227:         if (ik > i[j]) {
228:           SWAP2(i[k],i[j],Ii[k],Ii[j],tmp);
229:           ik = i[k];
230:         }
231:       }
232:     }
233:   } else {
234:     PetscSortIntWithArray_Private(i,Ii,n-1);
235:   }
236:   return(0);
237: }


240: #define SWAP3(a,b,c,d,e,f,t) {t=a;a=b;b=t;t=c;c=d;d=t;t=e;e=f;f=t;}

244: /*
245:    A simple version of quicksort; taken from Kernighan and Ritchie, page 87.
246:    Assumes 0 origin for v, number of elements = right+1 (right is index of
247:    right-most member).
248: */
249: static PetscErrorCode PetscSortIntWithArrayPair_Private(PetscInt *L,PetscInt *J, PetscInt *K, PetscInt right)
250: {
252:   PetscInt       i,vl,last,tmp;

255:   if (right <= 1) {
256:     if (right == 1) {
257:       if (L[0] > L[1]) SWAP3(L[0],L[1],J[0],J[1],K[0],K[1],tmp);
258:     }
259:     return(0);
260:   }
261:   SWAP3(L[0],L[right/2],J[0],J[right/2],K[0],K[right/2],tmp);
262:   vl   = L[0];
263:   last = 0;
264:   for (i=1; i<=right; i++) {
265:     if (L[i] < vl) {last++; SWAP3(L[last],L[i],J[last],J[i],K[last],K[i],tmp);}
266:   }
267:   SWAP3(L[0],L[last],J[0],J[last],K[0],K[last],tmp);
268:   PetscSortIntWithArrayPair_Private(L,J,K,last-1);
269:   PetscSortIntWithArrayPair_Private(L+last+1,J+last+1,K+last+1,right-(last+1));
270:   return(0);
271: }

275: /*@
276:    PetscSortIntWithArrayPair - Sorts an array of integers in place in increasing order;
277:        changes a pair of integer arrays to match the sorted first array.

279:    Not Collective

281:    Input Parameters:
282: +  n  - number of values
283: .  I  - array of integers
284: .  J  - second array of integers (first array of the pair)
285: -  K  - third array of integers  (second array of the pair)

287:    Level: intermediate

289:    Concepts: sorting^ints with array pair

291: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortIntWithArray()
292: @*/
293: PetscErrorCode  PetscSortIntWithArrayPair(PetscInt n,PetscInt *L,PetscInt *J, PetscInt *K)
294: {
296:   PetscInt       j,k,tmp,ik;

299:   if (n<8) {
300:     for (k=0; k<n; k++) {
301:       ik = L[k];
302:       for (j=k+1; j<n; j++) {
303:         if (ik > L[j]) {
304:           SWAP3(L[k],L[j],J[k],J[j],K[k],K[j],tmp);
305:           ik = L[k];
306:         }
307:       }
308:     }
309:   } else {
310:     PetscSortIntWithArrayPair_Private(L,J,K,n-1);
311:   }
312:   return(0);
313: }

317: /*
318:    A simple version of quicksort; taken from Kernighan and Ritchie, page 87.
319:    Assumes 0 origin for v, number of elements = right+1 (right is index of
320:    right-most member).
321: */
322: static void PetscSortMPIInt_Private(PetscMPIInt *v,PetscInt right)
323: {
324:   PetscInt          i,j;
325:   PetscMPIInt       pivot,tmp;

327:   if (right <= 1) {
328:     if (right == 1) {
329:       if (v[0] > v[1]) SWAP(v[0],v[1],tmp);
330:     }
331:     return;
332:   }
333:   i = MEDIAN(v,right);          /* Choose a pivot */
334:   SWAP(v[0],v[i],tmp);          /* Move it out of the way */
335:   pivot = v[0];
336:   for (i=0,j=right+1;; ) {
337:     while (++i < j && v[i] <= pivot) ; /* Scan from the left */
338:     while (v[--j] > pivot) ;           /* Scan from the right */
339:     if (i >= j) break;
340:     SWAP(v[i],v[j],tmp);
341:   }
342:   SWAP(v[0],v[j],tmp);          /* Put pivot back in place. */
343:   PetscSortMPIInt_Private(v,j-1);
344:   PetscSortMPIInt_Private(v+j+1,right-(j+1));
345: }

349: /*@
350:    PetscSortMPIInt - Sorts an array of MPI integers in place in increasing order.

352:    Not Collective

354:    Input Parameters:
355: +  n  - number of values
356: -  i  - array of integers

358:    Level: intermediate

360:    Concepts: sorting^ints

362: .seealso: PetscSortReal(), PetscSortIntWithPermutation()
363: @*/
364: PetscErrorCode  PetscSortMPIInt(PetscInt n,PetscMPIInt i[])
365: {
366:   PetscInt    j,k;
367:   PetscMPIInt tmp,ik;

370:   if (n<8) {
371:     for (k=0; k<n; k++) {
372:       ik = i[k];
373:       for (j=k+1; j<n; j++) {
374:         if (ik > i[j]) {
375:           SWAP(i[k],i[j],tmp);
376:           ik = i[k];
377:         }
378:       }
379:     }
380:   } else PetscSortMPIInt_Private(i,n-1);
381:   return(0);
382: }

386: /*@
387:    PetscSortRemoveDupsMPIInt - Sorts an array of MPI integers in place in increasing order removes all duplicate entries

389:    Not Collective

391:    Input Parameters:
392: +  n  - number of values
393: -  ii  - array of integers

395:    Output Parameter:
396: .  n - number of non-redundant values

398:    Level: intermediate

400:    Concepts: sorting^ints

402: .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt()
403: @*/
404: PetscErrorCode  PetscSortRemoveDupsMPIInt(PetscInt *n,PetscMPIInt ii[])
405: {
407:   PetscInt       i,s = 0,N = *n, b = 0;

410:   PetscSortMPIInt(N,ii);
411:   for (i=0; i<N-1; i++) {
412:     if (ii[b+s+1] != ii[b]) {
413:       ii[b+1] = ii[b+s+1]; b++;
414:     } else s++;
415:   }
416:   *n = N - s;
417:   return(0);
418: }

422: /*
423:    A simple version of quicksort; taken from Kernighan and Ritchie, page 87.
424:    Assumes 0 origin for v, number of elements = right+1 (right is index of
425:    right-most member).
426: */
427: static PetscErrorCode PetscSortMPIIntWithArray_Private(PetscMPIInt *v,PetscMPIInt *V,PetscMPIInt right)
428: {
430:   PetscMPIInt    i,vl,last,tmp;

433:   if (right <= 1) {
434:     if (right == 1) {
435:       if (v[0] > v[1]) SWAP2(v[0],v[1],V[0],V[1],tmp);
436:     }
437:     return(0);
438:   }
439:   SWAP2(v[0],v[right/2],V[0],V[right/2],tmp);
440:   vl   = v[0];
441:   last = 0;
442:   for (i=1; i<=right; i++) {
443:     if (v[i] < vl) {last++; SWAP2(v[last],v[i],V[last],V[i],tmp);}
444:   }
445:   SWAP2(v[0],v[last],V[0],V[last],tmp);
446:   PetscSortMPIIntWithArray_Private(v,V,last-1);
447:   PetscSortMPIIntWithArray_Private(v+last+1,V+last+1,right-(last+1));
448:   return(0);
449: }

453: /*@
454:    PetscSortMPIIntWithArray - Sorts an array of integers in place in increasing order;
455:        changes a second array to match the sorted first array.

457:    Not Collective

459:    Input Parameters:
460: +  n  - number of values
461: .  i  - array of integers
462: -  I - second array of integers

464:    Level: intermediate

466:    Concepts: sorting^ints with array

468: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt()
469: @*/
470: PetscErrorCode  PetscSortMPIIntWithArray(PetscMPIInt n,PetscMPIInt i[],PetscMPIInt Ii[])
471: {
473:   PetscMPIInt    j,k,tmp,ik;

476:   if (n<8) {
477:     for (k=0; k<n; k++) {
478:       ik = i[k];
479:       for (j=k+1; j<n; j++) {
480:         if (ik > i[j]) {
481:           SWAP2(i[k],i[j],Ii[k],Ii[j],tmp);
482:           ik = i[k];
483:         }
484:       }
485:     }
486:   } else {
487:     PetscSortMPIIntWithArray_Private(i,Ii,n-1);
488:   }
489:   return(0);
490: }

492: /* -----------------------------------------------------------------------*/
493: #define SWAP2IntScalar(a,b,c,d,t,ts) {t=a;a=b;b=t;ts=c;c=d;d=ts;}

497: /*
498:    Modified from PetscSortIntWithArray_Private().
499: */
500: static PetscErrorCode PetscSortIntWithScalarArray_Private(PetscInt *v,PetscScalar *V,PetscInt right)
501: {
503:   PetscInt       i,vl,last,tmp;
504:   PetscScalar    stmp;

507:   if (right <= 1) {
508:     if (right == 1) {
509:       if (v[0] > v[1]) SWAP2IntScalar(v[0],v[1],V[0],V[1],tmp,stmp);
510:     }
511:     return(0);
512:   }
513:   SWAP2IntScalar(v[0],v[right/2],V[0],V[right/2],tmp,stmp);
514:   vl   = v[0];
515:   last = 0;
516:   for (i=1; i<=right; i++) {
517:     if (v[i] < vl) {last++; SWAP2IntScalar(v[last],v[i],V[last],V[i],tmp,stmp);}
518:   }
519:   SWAP2IntScalar(v[0],v[last],V[0],V[last],tmp,stmp);
520:   PetscSortIntWithScalarArray_Private(v,V,last-1);
521:   PetscSortIntWithScalarArray_Private(v+last+1,V+last+1,right-(last+1));
522:   return(0);
523: }

527: /*@
528:    PetscSortIntWithScalarArray - Sorts an array of integers in place in increasing order;
529:        changes a second SCALAR array to match the sorted first INTEGER array.

531:    Not Collective

533:    Input Parameters:
534: +  n  - number of values
535: .  i  - array of integers
536: -  I - second array of scalars

538:    Level: intermediate

540:    Concepts: sorting^ints with array

542: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
543: @*/
544: PetscErrorCode  PetscSortIntWithScalarArray(PetscInt n,PetscInt i[],PetscScalar Ii[])
545: {
547:   PetscInt       j,k,tmp,ik;
548:   PetscScalar    stmp;

551:   if (n<8) {
552:     for (k=0; k<n; k++) {
553:       ik = i[k];
554:       for (j=k+1; j<n; j++) {
555:         if (ik > i[j]) {
556:           SWAP2IntScalar(i[k],i[j],Ii[k],Ii[j],tmp,stmp);
557:           ik = i[k];
558:         }
559:       }
560:     }
561:   } else {
562:     PetscSortIntWithScalarArray_Private(i,Ii,n-1);
563:   }
564:   return(0);
565: }

567: #define SWAP2IntData(a,b,c,d,t,td,siz)          \
568:   do {                                          \
569:   PetscErrorCode _ierr;                         \
570:   t=a;a=b;b=t;                                  \
571:   _PetscMemcpy(td,c,siz);CHKERRQ(_ierr); \
572:   _PetscMemcpy(c,d,siz);CHKERRQ(_ierr);  \
573:   _PetscMemcpy(d,td,siz);CHKERRQ(_ierr); \
574: } while(0)

578: /*
579:    Modified from PetscSortIntWithArray_Private().
580: */
581: static PetscErrorCode PetscSortIntWithDataArray_Private(PetscInt *v,char *V,PetscInt right,size_t size,void *work)
582: {
584:   PetscInt       i,vl,last,tmp;

587:   if (right <= 1) {
588:     if (right == 1) {
589:       if (v[0] > v[1]) SWAP2IntData(v[0],v[1],V,V+size,tmp,work,size);
590:     }
591:     return(0);
592:   }
593:   SWAP2IntData(v[0],v[right/2],V,V+size*(right/2),tmp,work,size);
594:   vl   = v[0];
595:   last = 0;
596:   for (i=1; i<=right; i++) {
597:     if (v[i] < vl) {last++; SWAP2IntData(v[last],v[i],V+size*last,V+size*i,tmp,work,size);}
598:   }
599:   SWAP2IntData(v[0],v[last],V,V + size*last,tmp,work,size);
600:   PetscSortIntWithDataArray_Private(v,V,last-1,size,work);
601:   PetscSortIntWithDataArray_Private(v+last+1,V+size*(last+1),right-(last+1),size,work);
602:   return(0);
603: }

607: /*@
608:    PetscSortIntWithDataArray - Sorts an array of integers in place in increasing order;
609:        changes a second array to match the sorted first INTEGER array.  Unlike other sort routines, the user must
610:        provide workspace (the size of an element in the data array) to use when sorting.

612:    Not Collective

614:    Input Parameters:
615: +  n  - number of values
616: .  i  - array of integers
617: .  Ii - second array of data
618: .  size - sizeof elements in the data array in bytes
619: -  work - workspace of "size" bytes used when sorting

621:    Level: intermediate

623:    Concepts: sorting^ints with array

625: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
626: @*/
627: PetscErrorCode  PetscSortIntWithDataArray(PetscInt n,PetscInt i[],void *Ii,size_t size,void *work)
628: {
629:   char           *V = (char *) Ii;
631:   PetscInt       j,k,tmp,ik;

634:   if (n<8) {
635:     for (k=0; k<n; k++) {
636:       ik = i[k];
637:       for (j=k+1; j<n; j++) {
638:         if (ik > i[j]) {
639:           SWAP2IntData(i[k],i[j],V+size*k,V+size*j,tmp,work,size);
640:           ik = i[k];
641:         }
642:       }
643:     }
644:   } else {
645:     PetscSortIntWithDataArray_Private(i,V,n-1,size,work);
646:   }
647:   return(0);
648: }


653: /*@
654:    PetscMergeIntArray -     Merges two SORTED integer arrays, removes duplicate elements.

656:    Not Collective

658:    Input Parameters:
659: +  an  - number of values in the first array
660: .  aI  - first sorted array of integers
661: .  bn  - number of values in the second array
662: -  bI  - second array of integers

664:    Output Parameters:
665: +  n   - number of values in the merged array
666: -  I   - merged sorted array, this is allocated if an array is not provided 

668:    Level: intermediate

670:    Concepts: merging^arrays

672: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
673: @*/
674: PetscErrorCode  PetscMergeIntArray(PetscInt an,const PetscInt *aI, PetscInt bn, const PetscInt *bI,  PetscInt *n, PetscInt **L)
675: {
677:   PetscInt       *L_ = *L, ak, bk, k;

679:   if (!L_) {
680:     PetscMalloc1(an+bn, L);
681:     L_   = *L;
682:   }
683:   k = ak = bk = 0;
684:   while (ak < an && bk < bn) {
685:     if (aI[ak] == bI[bk]) {
686:       L_[k] = aI[ak];
687:       ++ak;
688:       ++bk;
689:       ++k;
690:     } else if (aI[ak] < bI[bk]) {
691:       L_[k] = aI[ak];
692:       ++ak;
693:       ++k;
694:     } else {
695:       L_[k] = bI[bk];
696:       ++bk;
697:       ++k;
698:     }
699:   }
700:   if (ak < an) {
701:     PetscMemcpy(L_+k,aI+ak,(an-ak)*sizeof(PetscInt));
702:     k   += (an-ak);
703:   }
704:   if (bk < bn) {
705:     PetscMemcpy(L_+k,bI+bk,(bn-bk)*sizeof(PetscInt));
706:     k   += (bn-bk);
707:   }
708:   *n = k;
709:   return(0);
710: }

714: /*@
715:    PetscMergeIntArrayPair -     Merges two SORTED integer arrays that share NO common values along with an additional array of integers.
716:                                 The additional arrays are the same length as sorted arrays and are merged
717:                                 in the order determined by the merging of the sorted pair.

719:    Not Collective

721:    Input Parameters:
722: +  an  - number of values in the first array
723: .  aI  - first sorted array of integers
724: .  aJ  - first additional array of integers
725: .  bn  - number of values in the second array
726: .  bI  - second array of integers
727: -  bJ  - second additional of integers

729:    Output Parameters:
730: +  n   - number of values in the merged array (== an + bn)
731: .  I   - merged sorted array
732: -  J   - merged additional array

734:    Level: intermediate

736:    Concepts: merging^arrays

738: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
739: @*/
740: PetscErrorCode  PetscMergeIntArrayPair(PetscInt an,const PetscInt *aI, const PetscInt *aJ, PetscInt bn, const PetscInt *bI, const PetscInt *bJ, PetscInt *n, PetscInt **L, PetscInt **J)
741: {
743:   PetscInt       n_, *L_ = *L, *J_= *J, ak, bk, k;

745:   n_ = an + bn;
746:   *n = n_;
747:   if (!L_) {
748:     PetscMalloc1(n_, L);
749:     L_   = *L;
750:   }
751:   if (!J_) {
752:     PetscMalloc1(n_, &J_);
753:     J_   = *J;
754:   }
755:   k = ak = bk = 0;
756:   while (ak < an && bk < bn) {
757:     if (aI[ak] <= bI[bk]) {
758:       L_[k] = aI[ak];
759:       J_[k] = aJ[ak];
760:       ++ak;
761:       ++k;
762:     } else {
763:       L_[k] = bI[bk];
764:       J_[k] = bJ[bk];
765:       ++bk;
766:       ++k;
767:     }
768:   }
769:   if (ak < an) {
770:     PetscMemcpy(L_+k,aI+ak,(an-ak)*sizeof(PetscInt));
771:     PetscMemcpy(J_+k,aJ+ak,(an-ak)*sizeof(PetscInt));
772:     k   += (an-ak);
773:   }
774:   if (bk < bn) {
775:     PetscMemcpy(L_+k,bI+bk,(bn-bk)*sizeof(PetscInt));
776:     PetscMemcpy(J_+k,bJ+bk,(bn-bk)*sizeof(PetscInt));
777:   }
778:   return(0);
779: }


784: /*@
785:    PetscProcessTree - Prepares tree data to be displayed graphically

787:    Not Collective

789:    Input Parameters:
790: +  n  - number of values
791: .  mask - indicates those entries in the tree, location 0 is always masked
792: -  parentid - indicates the parent of each entry

794:    Output Parameters:
795: +  Nlevels - the number of levels
796: .  Level - for each node tells its level
797: .  Levelcnts - the number of nodes on each level
798: .  Idbylevel - a list of ids on each of the levels, first level followed by second etc
799: -  Column - for each id tells its column index

801:    Level: intermediate


804: .seealso: PetscSortReal(), PetscSortIntWithPermutation()
805: @*/
806: PetscErrorCode  PetscProcessTree(PetscInt n,const PetscBool mask[],const PetscInt parentid[],PetscInt *Nlevels,PetscInt **Level,PetscInt **Levelcnt,PetscInt **Idbylevel,PetscInt **Column)
807: {
808:   PetscInt       i,j,cnt,nmask = 0,nlevels = 0,*level,*levelcnt,levelmax = 0,*workid,*workparentid,tcnt = 0,*idbylevel,*column;
810:   PetscBool      done = PETSC_FALSE;

813:   if (!mask[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mask of 0th location must be set");
814:   for (i=0; i<n; i++) {
815:     if (mask[i]) continue;
816:     if (parentid[i]  == i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Node labeled as own parent");
817:     if (parentid[i] && mask[parentid[i]]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Parent is masked");
818:   }

820:   for (i=0; i<n; i++) {
821:     if (!mask[i]) nmask++;
822:   }

824:   /* determine the level in the tree of each node */
825:   PetscCalloc1(n,&level);

827:   level[0] = 1;
828:   while (!done) {
829:     done = PETSC_TRUE;
830:     for (i=0; i<n; i++) {
831:       if (mask[i]) continue;
832:       if (!level[i] && level[parentid[i]]) level[i] = level[parentid[i]] + 1;
833:       else if (!level[i]) done = PETSC_FALSE;
834:     }
835:   }
836:   for (i=0; i<n; i++) {
837:     level[i]--;
838:     nlevels = PetscMax(nlevels,level[i]);
839:   }

841:   /* count the number of nodes on each level and its max */
842:   PetscCalloc1(nlevels,&levelcnt);
843:   for (i=0; i<n; i++) {
844:     if (mask[i]) continue;
845:     levelcnt[level[i]-1]++;
846:   }
847:   for (i=0; i<nlevels;i++) levelmax = PetscMax(levelmax,levelcnt[i]);

849:   /* for each level sort the ids by the parent id */
850:   PetscMalloc2(levelmax,&workid,levelmax,&workparentid);
851:   PetscMalloc1(nmask,&idbylevel);
852:   for (j=1; j<=nlevels;j++) {
853:     cnt = 0;
854:     for (i=0; i<n; i++) {
855:       if (mask[i]) continue;
856:       if (level[i] != j) continue;
857:       workid[cnt]         = i;
858:       workparentid[cnt++] = parentid[i];
859:     }
860:     /*  PetscIntView(cnt,workparentid,0);
861:     PetscIntView(cnt,workid,0);
862:     PetscSortIntWithArray(cnt,workparentid,workid);
863:     PetscIntView(cnt,workparentid,0);
864:     PetscIntView(cnt,workid,0);*/
865:     PetscMemcpy(idbylevel+tcnt,workid,cnt*sizeof(PetscInt));
866:     tcnt += cnt;
867:   }
868:   if (tcnt != nmask) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent count of unmasked nodes");
869:   PetscFree2(workid,workparentid);

871:   /* for each node list its column */
872:   PetscMalloc1(n,&column);
873:   cnt = 0;
874:   for (j=0; j<nlevels; j++) {
875:     for (i=0; i<levelcnt[j]; i++) {
876:       column[idbylevel[cnt++]] = i;
877:     }
878:   }

880:   *Nlevels   = nlevels;
881:   *Level     = level;
882:   *Levelcnt  = levelcnt;
883:   *Idbylevel = idbylevel;
884:   *Column    = column;
885:   return(0);
886: }