Actual source code: sorti.c

petsc-3.5.1 2014-07-24
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: }


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

573:    Not Collective

575:    Input Parameters:
576: +  an  - number of values in the first array
577: .  aI  - first sorted array of integers
578: .  bn  - number of values in the second array
579: -  bI  - second array of integers

581:    Output Parameters:
582: +  n   - number of values in the merged array
583: -  I   - merged sorted array, this is allocated if an array is not provided 

585:    Level: intermediate

587:    Concepts: merging^arrays

589: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
590: @*/
591: PetscErrorCode  PetscMergeIntArray(PetscInt an,const PetscInt *aI, PetscInt bn, const PetscInt *bI,  PetscInt *n, PetscInt **L)
592: {
594:   PetscInt       *L_ = *L, ak, bk, k;

596:   if (!L_) {
597:     PetscMalloc1(an+bn, L);
598:     L_   = *L;
599:   }
600:   k = ak = bk = 0;
601:   while (ak < an && bk < bn) {
602:     if (aI[ak] == bI[bk]) {
603:       L_[k] = aI[ak];
604:       ++ak;
605:       ++bk;
606:       ++k;
607:     } else if (aI[ak] < bI[bk]) {
608:       L_[k] = aI[ak];
609:       ++ak;
610:       ++k;
611:     } else {
612:       L_[k] = bI[bk];
613:       ++bk;
614:       ++k;
615:     }
616:   }
617:   if (ak < an) {
618:     PetscMemcpy(L_+k,aI+ak,(an-ak)*sizeof(PetscInt));
619:     k   += (an-ak);
620:   }
621:   if (bk < bn) {
622:     PetscMemcpy(L_+k,bI+bk,(bn-bk)*sizeof(PetscInt));
623:     k   += (bn-bk);
624:   }
625:   *n = k;
626:   return(0);
627: }

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

636:    Not Collective

638:    Input Parameters:
639: +  an  - number of values in the first array
640: .  aI  - first sorted array of integers
641: .  aJ  - first additional array of integers
642: .  bn  - number of values in the second array
643: .  bI  - second array of integers
644: -  bJ  - second additional of integers

646:    Output Parameters:
647: +  n   - number of values in the merged array (== an + bn)
648: .  I   - merged sorted array
649: -  J   - merged additional array

651:    Level: intermediate

653:    Concepts: merging^arrays

655: .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray()
656: @*/
657: PetscErrorCode  PetscMergeIntArrayPair(PetscInt an,const PetscInt *aI, const PetscInt *aJ, PetscInt bn, const PetscInt *bI, const PetscInt *bJ, PetscInt *n, PetscInt **L, PetscInt **J)
658: {
660:   PetscInt       n_, *L_ = *L, *J_= *J, ak, bk, k;

662:   n_ = an + bn;
663:   *n = n_;
664:   if (!L_) {
665:     PetscMalloc1(n_, L);
666:     L_   = *L;
667:   }
668:   if (!J_) {
669:     PetscMalloc1(n_, &J_);
670:     J_   = *J;
671:   }
672:   k = ak = bk = 0;
673:   while (ak < an && bk < bn) {
674:     if (aI[ak] <= bI[bk]) {
675:       L_[k] = aI[ak];
676:       J_[k] = aJ[ak];
677:       ++ak;
678:       ++k;
679:     } else {
680:       L_[k] = bI[bk];
681:       J_[k] = bJ[bk];
682:       ++bk;
683:       ++k;
684:     }
685:   }
686:   if (ak < an) {
687:     PetscMemcpy(L_+k,aI+ak,(an-ak)*sizeof(PetscInt));
688:     PetscMemcpy(J_+k,aJ+ak,(an-ak)*sizeof(PetscInt));
689:     k   += (an-ak);
690:   }
691:   if (bk < bn) {
692:     PetscMemcpy(L_+k,bI+bk,(bn-bk)*sizeof(PetscInt));
693:     PetscMemcpy(J_+k,bJ+bk,(bn-bk)*sizeof(PetscInt));
694:   }
695:   return(0);
696: }


701: /*@
702:    PetscProcessTree - Prepares tree data to be displayed graphically

704:    Not Collective

706:    Input Parameters:
707: +  n  - number of values
708: .  mask - indicates those entries in the tree, location 0 is always masked
709: -  parentid - indicates the parent of each entry

711:    Output Parameters:
712: +  Nlevels - the number of levels
713: .  Level - for each node tells its level
714: .  Levelcnts - the number of nodes on each level
715: .  Idbylevel - a list of ids on each of the levels, first level followed by second etc
716: -  Column - for each id tells its column index

718:    Level: intermediate


721: .seealso: PetscSortReal(), PetscSortIntWithPermutation()
722: @*/
723: PetscErrorCode  PetscProcessTree(PetscInt n,const PetscBool mask[],const PetscInt parentid[],PetscInt *Nlevels,PetscInt **Level,PetscInt **Levelcnt,PetscInt **Idbylevel,PetscInt **Column)
724: {
725:   PetscInt       i,j,cnt,nmask = 0,nlevels = 0,*level,*levelcnt,levelmax = 0,*workid,*workparentid,tcnt = 0,*idbylevel,*column;
727:   PetscBool      done = PETSC_FALSE;

730:   if (!mask[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mask of 0th location must be set");
731:   for (i=0; i<n; i++) {
732:     if (mask[i]) continue;
733:     if (parentid[i]  == i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Node labeled as own parent");
734:     if (parentid[i] && mask[parentid[i]]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Parent is masked");
735:   }

737:   for (i=0; i<n; i++) {
738:     if (!mask[i]) nmask++;
739:   }

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

744:   level[0] = 1;
745:   while (!done) {
746:     done = PETSC_TRUE;
747:     for (i=0; i<n; i++) {
748:       if (mask[i]) continue;
749:       if (!level[i] && level[parentid[i]]) level[i] = level[parentid[i]] + 1;
750:       else if (!level[i]) done = PETSC_FALSE;
751:     }
752:   }
753:   for (i=0; i<n; i++) {
754:     level[i]--;
755:     nlevels = PetscMax(nlevels,level[i]);
756:   }

758:   /* count the number of nodes on each level and its max */
759:   PetscCalloc1(nlevels,&levelcnt);
760:   for (i=0; i<n; i++) {
761:     if (mask[i]) continue;
762:     levelcnt[level[i]-1]++;
763:   }
764:   for (i=0; i<nlevels;i++) levelmax = PetscMax(levelmax,levelcnt[i]);

766:   /* for each level sort the ids by the parent id */
767:   PetscMalloc2(levelmax,&workid,levelmax,&workparentid);
768:   PetscMalloc1(nmask,&idbylevel);
769:   for (j=1; j<=nlevels;j++) {
770:     cnt = 0;
771:     for (i=0; i<n; i++) {
772:       if (mask[i]) continue;
773:       if (level[i] != j) continue;
774:       workid[cnt]         = i;
775:       workparentid[cnt++] = parentid[i];
776:     }
777:     /*  PetscIntView(cnt,workparentid,0);
778:     PetscIntView(cnt,workid,0);
779:     PetscSortIntWithArray(cnt,workparentid,workid);
780:     PetscIntView(cnt,workparentid,0);
781:     PetscIntView(cnt,workid,0);*/
782:     PetscMemcpy(idbylevel+tcnt,workid,cnt*sizeof(PetscInt));
783:     tcnt += cnt;
784:   }
785:   if (tcnt != nmask) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent count of unmasked nodes");
786:   PetscFree2(workid,workparentid);

788:   /* for each node list its column */
789:   PetscMalloc1(n,&column);
790:   cnt = 0;
791:   for (j=0; j<nlevels; j++) {
792:     for (i=0; i<levelcnt[j]; i++) {
793:       column[idbylevel[cnt++]] = i;
794:     }
795:   }

797:   *Nlevels   = nlevels;
798:   *Level     = level;
799:   *Levelcnt  = levelcnt;
800:   *Idbylevel = idbylevel;
801:   *Column    = column;
802:   return(0);
803: }