1 2 /* 3 This file contains routines for sorting integers. Values are sorted in place. 4 */ 5 #include <petscsys.h> /*I "petscsys.h" I*/ 6 7 #define SWAP(a,b,t) {t=a;a=b;b=t;} 8 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))) 17 18 #define MEDIAN(v,right) MEDIAN3(v,right/4,right/2,right/4*3) 19 20 /* -----------------------------------------------------------------------*/ 21 22 #undef __FUNCT__ 23 #define __FUNCT__ "PetscSortInt_Private" 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; 32 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 } 52 53 #undef __FUNCT__ 54 #define __FUNCT__ "PetscSortInt" 55 /*@ 56 PetscSortInt - Sorts an array of integers in place in increasing order. 57 58 Not Collective 59 60 Input Parameters: 61 + n - number of values 62 - i - array of integers 63 64 Level: intermediate 65 66 Concepts: sorting^ints 67 68 .seealso: PetscSortReal(), PetscSortIntWithPermutation() 69 @*/ 70 PetscErrorCode PetscSortInt(PetscInt n,PetscInt i[]) 71 { 72 PetscInt j,k,tmp,ik; 73 74 PetscFunctionBegin; 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 { 86 PetscSortInt_Private(i,n-1); 87 } 88 PetscFunctionReturn(0); 89 } 90 91 #undef __FUNCT__ 92 #define __FUNCT__ "PetscSortRemoveDupsInt" 93 /*@ 94 PetscSortRemoveDupsInt - Sorts an array of integers in place in increasing order removes all duplicate entries 95 96 Not Collective 97 98 Input Parameters: 99 + n - number of values 100 - ii - array of integers 101 102 Output Parameter: 103 . n - number of non-redundant values 104 105 Level: intermediate 106 107 Concepts: sorting^ints 108 109 .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt() 110 @*/ 111 PetscErrorCode PetscSortRemoveDupsInt(PetscInt *n,PetscInt ii[]) 112 { 113 PetscErrorCode ierr; 114 PetscInt i,s = 0,N = *n, b = 0; 115 116 PetscFunctionBegin; 117 ierr = PetscSortInt(N,ii);CHKERRQ(ierr); 118 for (i=0; i<N-1; i++) { 119 if (ii[b+s+1] != ii[b]) {ii[b+1] = ii[b+s+1]; b++;} 120 else s++; 121 } 122 *n = N - s; 123 PetscFunctionReturn(0); 124 } 125 126 #undef __FUNCT__ 127 #define __FUNCT__ "PetscFindInt" 128 /*@ 129 PetscFindInt - Finds integer in a sorted array of integers 130 131 Not Collective 132 133 Input Parameters: 134 + key - the integer to locate 135 . n - number of values in the array 136 - ii - array of integers 137 138 Output Parameter: 139 . loc - the location if found, otherwise -(slot+1) where slot is the place the value would go 140 141 Level: intermediate 142 143 Concepts: sorting^ints 144 145 .seealso: PetscSortInt(), PetscSortIntWithArray(), PetscSortRemoveDupsInt() 146 @*/ 147 PetscErrorCode PetscFindInt(PetscInt key, PetscInt n, const PetscInt ii[], PetscInt *loc) 148 { 149 PetscInt lo = 0,hi = n; 150 151 PetscFunctionBegin; 152 PetscValidPointer(loc,4); 153 if (!n) {*loc = -1; PetscFunctionReturn(0);} 154 PetscValidPointer(ii,3); 155 while (hi - lo > 1) { 156 PetscInt mid = lo + (hi - lo)/2; 157 if (key < ii[mid]) hi = mid; 158 else lo = mid; 159 } 160 *loc = key == ii[lo] ? lo : -(lo + (key > ii[lo]) + 1); 161 PetscFunctionReturn(0); 162 } 163 164 165 /* -----------------------------------------------------------------------*/ 166 #define SWAP2(a,b,c,d,t) {t=a;a=b;b=t;t=c;c=d;d=t;} 167 168 #undef __FUNCT__ 169 #define __FUNCT__ "PetscSortIntWithArray_Private" 170 /* 171 A simple version of quicksort; taken from Kernighan and Ritchie, page 87. 172 Assumes 0 origin for v, number of elements = right+1 (right is index of 173 right-most member). 174 */ 175 static PetscErrorCode PetscSortIntWithArray_Private(PetscInt *v,PetscInt *V,PetscInt right) 176 { 177 PetscErrorCode ierr; 178 PetscInt i,vl,last,tmp; 179 180 PetscFunctionBegin; 181 if (right <= 1) { 182 if (right == 1) { 183 if (v[0] > v[1]) SWAP2(v[0],v[1],V[0],V[1],tmp); 184 } 185 PetscFunctionReturn(0); 186 } 187 SWAP2(v[0],v[right/2],V[0],V[right/2],tmp); 188 vl = v[0]; 189 last = 0; 190 for (i=1; i<=right; i++) { 191 if (v[i] < vl) {last++; SWAP2(v[last],v[i],V[last],V[i],tmp);} 192 } 193 SWAP2(v[0],v[last],V[0],V[last],tmp); 194 ierr = PetscSortIntWithArray_Private(v,V,last-1);CHKERRQ(ierr); 195 ierr = PetscSortIntWithArray_Private(v+last+1,V+last+1,right-(last+1));CHKERRQ(ierr); 196 PetscFunctionReturn(0); 197 } 198 199 #undef __FUNCT__ 200 #define __FUNCT__ "PetscSortIntWithArray" 201 /*@ 202 PetscSortIntWithArray - Sorts an array of integers in place in increasing order; 203 changes a second array to match the sorted first array. 204 205 Not Collective 206 207 Input Parameters: 208 + n - number of values 209 . i - array of integers 210 - I - second array of integers 211 212 Level: intermediate 213 214 Concepts: sorting^ints with array 215 216 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt() 217 @*/ 218 PetscErrorCode PetscSortIntWithArray(PetscInt n,PetscInt i[],PetscInt Ii[]) 219 { 220 PetscErrorCode ierr; 221 PetscInt j,k,tmp,ik; 222 223 PetscFunctionBegin; 224 if (n<8) { 225 for (k=0; k<n; k++) { 226 ik = i[k]; 227 for (j=k+1; j<n; j++) { 228 if (ik > i[j]) { 229 SWAP2(i[k],i[j],Ii[k],Ii[j],tmp); 230 ik = i[k]; 231 } 232 } 233 } 234 } else { 235 ierr = PetscSortIntWithArray_Private(i,Ii,n-1);CHKERRQ(ierr); 236 } 237 PetscFunctionReturn(0); 238 } 239 240 241 #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;} 242 243 #undef __FUNCT__ 244 #define __FUNCT__ "PetscSortIntWithArrayPair_Private" 245 /* 246 A simple version of quicksort; taken from Kernighan and Ritchie, page 87. 247 Assumes 0 origin for v, number of elements = right+1 (right is index of 248 right-most member). 249 */ 250 static PetscErrorCode PetscSortIntWithArrayPair_Private(PetscInt *L,PetscInt *J, PetscInt *K, PetscInt right) 251 { 252 PetscErrorCode ierr; 253 PetscInt i,vl,last,tmp; 254 255 PetscFunctionBegin; 256 if (right <= 1) { 257 if (right == 1) { 258 if (L[0] > L[1]) SWAP3(L[0],L[1],J[0],J[1],K[0],K[1],tmp); 259 } 260 PetscFunctionReturn(0); 261 } 262 SWAP3(L[0],L[right/2],J[0],J[right/2],K[0],K[right/2],tmp); 263 vl = L[0]; 264 last = 0; 265 for (i=1; i<=right; i++) { 266 if (L[i] < vl) {last++; SWAP3(L[last],L[i],J[last],J[i],K[last],K[i],tmp);} 267 } 268 SWAP3(L[0],L[last],J[0],J[last],K[0],K[last],tmp); 269 ierr = PetscSortIntWithArrayPair_Private(L,J,K,last-1);CHKERRQ(ierr); 270 ierr = PetscSortIntWithArrayPair_Private(L+last+1,J+last+1,K+last+1,right-(last+1));CHKERRQ(ierr); 271 PetscFunctionReturn(0); 272 } 273 274 #undef __FUNCT__ 275 #define __FUNCT__ "PetscSortIntWithArrayPair" 276 /*@ 277 PetscSortIntWithArrayPair - Sorts an array of integers in place in increasing order; 278 changes a pair of integer arrays to match the sorted first array. 279 280 Not Collective 281 282 Input Parameters: 283 + n - number of values 284 . I - array of integers 285 . J - second array of integers (first array of the pair) 286 - K - third array of integers (second array of the pair) 287 288 Level: intermediate 289 290 Concepts: sorting^ints with array pair 291 292 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortIntWithArray() 293 @*/ 294 PetscErrorCode PetscSortIntWithArrayPair(PetscInt n,PetscInt *L,PetscInt *J, PetscInt *K) 295 { 296 PetscErrorCode ierr; 297 PetscInt j,k,tmp,ik; 298 299 PetscFunctionBegin; 300 if (n<8) { 301 for (k=0; k<n; k++) { 302 ik = L[k]; 303 for (j=k+1; j<n; j++) { 304 if (ik > L[j]) { 305 SWAP3(L[k],L[j],J[k],J[j],K[k],K[j],tmp); 306 ik = L[k]; 307 } 308 } 309 } 310 } else { 311 ierr = PetscSortIntWithArrayPair_Private(L,J,K,n-1);CHKERRQ(ierr); 312 } 313 PetscFunctionReturn(0); 314 } 315 316 #undef __FUNCT__ 317 #define __FUNCT__ "PetscSortMPIInt_Private" 318 /* 319 A simple version of quicksort; taken from Kernighan and Ritchie, page 87. 320 Assumes 0 origin for v, number of elements = right+1 (right is index of 321 right-most member). 322 */ 323 static void PetscSortMPIInt_Private(PetscMPIInt *v,PetscInt right) 324 { 325 PetscInt i,j; 326 PetscMPIInt pivot,tmp; 327 328 if (right <= 1) { 329 if (right == 1) { 330 if (v[0] > v[1]) SWAP(v[0],v[1],tmp); 331 } 332 return; 333 } 334 i = MEDIAN(v,right); /* Choose a pivot */ 335 SWAP(v[0],v[i],tmp); /* Move it out of the way */ 336 pivot = v[0]; 337 for (i=0,j=right+1;;) { 338 while (++i < j && v[i] <= pivot); /* Scan from the left */ 339 while (v[--j] > pivot); /* Scan from the right */ 340 if (i >= j) break; 341 SWAP(v[i],v[j],tmp); 342 } 343 SWAP(v[0],v[j],tmp); /* Put pivot back in place. */ 344 PetscSortMPIInt_Private(v,j-1); 345 PetscSortMPIInt_Private(v+j+1,right-(j+1)); 346 } 347 348 #undef __FUNCT__ 349 #define __FUNCT__ "PetscSortMPIInt" 350 /*@ 351 PetscSortMPIInt - Sorts an array of MPI integers in place in increasing order. 352 353 Not Collective 354 355 Input Parameters: 356 + n - number of values 357 - i - array of integers 358 359 Level: intermediate 360 361 Concepts: sorting^ints 362 363 .seealso: PetscSortReal(), PetscSortIntWithPermutation() 364 @*/ 365 PetscErrorCode PetscSortMPIInt(PetscInt n,PetscMPIInt i[]) 366 { 367 PetscInt j,k; 368 PetscMPIInt tmp,ik; 369 370 PetscFunctionBegin; 371 if (n<8) { 372 for (k=0; k<n; k++) { 373 ik = i[k]; 374 for (j=k+1; j<n; j++) { 375 if (ik > i[j]) { 376 SWAP(i[k],i[j],tmp); 377 ik = i[k]; 378 } 379 } 380 } 381 } else { 382 PetscSortMPIInt_Private(i,n-1); 383 } 384 PetscFunctionReturn(0); 385 } 386 387 #undef __FUNCT__ 388 #define __FUNCT__ "PetscSortRemoveDupsMPIInt" 389 /*@ 390 PetscSortRemoveDupsMPIInt - Sorts an array of MPI integers in place in increasing order removes all duplicate entries 391 392 Not Collective 393 394 Input Parameters: 395 + n - number of values 396 - ii - array of integers 397 398 Output Parameter: 399 . n - number of non-redundant values 400 401 Level: intermediate 402 403 Concepts: sorting^ints 404 405 .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt() 406 @*/ 407 PetscErrorCode PetscSortRemoveDupsMPIInt(PetscInt *n,PetscMPIInt ii[]) 408 { 409 PetscErrorCode ierr; 410 PetscInt i,s = 0,N = *n, b = 0; 411 412 PetscFunctionBegin; 413 ierr = PetscSortMPIInt(N,ii);CHKERRQ(ierr); 414 for (i=0; i<N-1; i++) { 415 if (ii[b+s+1] != ii[b]) {ii[b+1] = ii[b+s+1]; b++;} 416 else s++; 417 } 418 *n = N - s; 419 PetscFunctionReturn(0); 420 } 421 422 #undef __FUNCT__ 423 #define __FUNCT__ "PetscSortMPIIntWithArray_Private" 424 /* 425 A simple version of quicksort; taken from Kernighan and Ritchie, page 87. 426 Assumes 0 origin for v, number of elements = right+1 (right is index of 427 right-most member). 428 */ 429 static PetscErrorCode PetscSortMPIIntWithArray_Private(PetscMPIInt *v,PetscMPIInt *V,PetscMPIInt right) 430 { 431 PetscErrorCode ierr; 432 PetscMPIInt i,vl,last,tmp; 433 434 PetscFunctionBegin; 435 if (right <= 1) { 436 if (right == 1) { 437 if (v[0] > v[1]) SWAP2(v[0],v[1],V[0],V[1],tmp); 438 } 439 PetscFunctionReturn(0); 440 } 441 SWAP2(v[0],v[right/2],V[0],V[right/2],tmp); 442 vl = v[0]; 443 last = 0; 444 for (i=1; i<=right; i++) { 445 if (v[i] < vl) {last++; SWAP2(v[last],v[i],V[last],V[i],tmp);} 446 } 447 SWAP2(v[0],v[last],V[0],V[last],tmp); 448 ierr = PetscSortMPIIntWithArray_Private(v,V,last-1);CHKERRQ(ierr); 449 ierr = PetscSortMPIIntWithArray_Private(v+last+1,V+last+1,right-(last+1));CHKERRQ(ierr); 450 PetscFunctionReturn(0); 451 } 452 453 #undef __FUNCT__ 454 #define __FUNCT__ "PetscSortMPIIntWithArray" 455 /*@ 456 PetscSortMPIIntWithArray - Sorts an array of integers in place in increasing order; 457 changes a second array to match the sorted first array. 458 459 Not Collective 460 461 Input Parameters: 462 + n - number of values 463 . i - array of integers 464 - I - second array of integers 465 466 Level: intermediate 467 468 Concepts: sorting^ints with array 469 470 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt() 471 @*/ 472 PetscErrorCode PetscSortMPIIntWithArray(PetscMPIInt n,PetscMPIInt i[],PetscMPIInt Ii[]) 473 { 474 PetscErrorCode ierr; 475 PetscMPIInt j,k,tmp,ik; 476 477 PetscFunctionBegin; 478 if (n<8) { 479 for (k=0; k<n; k++) { 480 ik = i[k]; 481 for (j=k+1; j<n; j++) { 482 if (ik > i[j]) { 483 SWAP2(i[k],i[j],Ii[k],Ii[j],tmp); 484 ik = i[k]; 485 } 486 } 487 } 488 } else { 489 ierr = PetscSortMPIIntWithArray_Private(i,Ii,n-1);CHKERRQ(ierr); 490 } 491 PetscFunctionReturn(0); 492 } 493 494 /* -----------------------------------------------------------------------*/ 495 #define SWAP2IntScalar(a,b,c,d,t,ts) {t=a;a=b;b=t;ts=c;c=d;d=ts;} 496 497 #undef __FUNCT__ 498 #define __FUNCT__ "PetscSortIntWithScalarArray_Private" 499 /* 500 Modified from PetscSortIntWithArray_Private(). 501 */ 502 static PetscErrorCode PetscSortIntWithScalarArray_Private(PetscInt *v,PetscScalar *V,PetscInt right) 503 { 504 PetscErrorCode ierr; 505 PetscInt i,vl,last,tmp; 506 PetscScalar stmp; 507 508 PetscFunctionBegin; 509 if (right <= 1) { 510 if (right == 1) { 511 if (v[0] > v[1]) SWAP2IntScalar(v[0],v[1],V[0],V[1],tmp,stmp); 512 } 513 PetscFunctionReturn(0); 514 } 515 SWAP2IntScalar(v[0],v[right/2],V[0],V[right/2],tmp,stmp); 516 vl = v[0]; 517 last = 0; 518 for (i=1; i<=right; i++) { 519 if (v[i] < vl) {last++; SWAP2IntScalar(v[last],v[i],V[last],V[i],tmp,stmp);} 520 } 521 SWAP2IntScalar(v[0],v[last],V[0],V[last],tmp,stmp); 522 ierr = PetscSortIntWithScalarArray_Private(v,V,last-1);CHKERRQ(ierr); 523 ierr = PetscSortIntWithScalarArray_Private(v+last+1,V+last+1,right-(last+1));CHKERRQ(ierr); 524 PetscFunctionReturn(0); 525 } 526 527 #undef __FUNCT__ 528 #define __FUNCT__ "PetscSortIntWithScalarArray" 529 /*@ 530 PetscSortIntWithScalarArray - Sorts an array of integers in place in increasing order; 531 changes a second SCALAR array to match the sorted first INTEGER array. 532 533 Not Collective 534 535 Input Parameters: 536 + n - number of values 537 . i - array of integers 538 - I - second array of scalars 539 540 Level: intermediate 541 542 Concepts: sorting^ints with array 543 544 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray() 545 @*/ 546 PetscErrorCode PetscSortIntWithScalarArray(PetscInt n,PetscInt i[],PetscScalar Ii[]) 547 { 548 PetscErrorCode ierr; 549 PetscInt j,k,tmp,ik; 550 PetscScalar stmp; 551 552 PetscFunctionBegin; 553 if (n<8) { 554 for (k=0; k<n; k++) { 555 ik = i[k]; 556 for (j=k+1; j<n; j++) { 557 if (ik > i[j]) { 558 SWAP2IntScalar(i[k],i[j],Ii[k],Ii[j],tmp,stmp); 559 ik = i[k]; 560 } 561 } 562 } 563 } else { 564 ierr = PetscSortIntWithScalarArray_Private(i,Ii,n-1);CHKERRQ(ierr); 565 } 566 PetscFunctionReturn(0); 567 } 568 569 570 571 #undef __FUNCT__ 572 #define __FUNCT__ "PetscMergeIntArrayPair" 573 /*@ 574 PetscMergeIntArrayPair - Merges two SORTED integer arrays along with an additional array of integers. 575 The additional arrays are the same length as sorted arrays and are merged 576 in the order determined by the merging of the sorted pair. 577 578 Not Collective 579 580 Input Parameters: 581 + an - number of values in the first array 582 . aI - first sorted array of integers 583 . aJ - first additional array of integers 584 . bn - number of values in the second array 585 . bI - second array of integers 586 - bJ - second additional of integers 587 588 Output Parameters: 589 + n - number of values in the merged array (== an + bn) 590 . I - merged sorted array 591 - J - merged additional array 592 593 Level: intermediate 594 595 Concepts: merging^arrays 596 597 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray() 598 @*/ 599 PetscErrorCode PetscMergeIntArrayPair(PetscInt an,const PetscInt *aI, const PetscInt *aJ, PetscInt bn, const PetscInt *bI, const PetscInt *bJ, PetscInt *n, PetscInt **L, PetscInt **J) 600 { 601 PetscErrorCode ierr; 602 PetscInt n_, *L_ = *L, *J_= *J, ak, bk, k; 603 604 n_ = an + bn; 605 *n = n_; 606 if (!L_) { 607 ierr = PetscMalloc(n_*sizeof(PetscInt), L);CHKERRQ(ierr); 608 L_ = *L; 609 } 610 if (!J_){ 611 ierr = PetscMalloc(n_*sizeof(PetscInt), &J_);CHKERRQ(ierr); 612 J_ = *J; 613 } 614 k = ak = bk = 0; 615 while(ak < an && bk < bn) { 616 if (aI[ak] <= bI[bk]) { 617 L_[k] = aI[ak]; 618 J_[k] = aJ[ak]; 619 ++ak; 620 ++k; 621 } 622 else { 623 L_[k] = bI[bk]; 624 J_[k] = bJ[bk]; 625 ++bk; 626 ++k; 627 } 628 } 629 if (ak < an) { 630 ierr = PetscMemcpy(L_+k,aI+ak,(an-ak)*sizeof(PetscInt));CHKERRQ(ierr); 631 ierr = PetscMemcpy(J_+k,aJ+ak,(an-ak)*sizeof(PetscInt));CHKERRQ(ierr); 632 k += (an-ak); 633 } 634 if (bk < bn) { 635 ierr = PetscMemcpy(L_+k,bI+bk,(bn-bk)*sizeof(PetscInt));CHKERRQ(ierr); 636 ierr = PetscMemcpy(J_+k,bJ+bk,(bn-bk)*sizeof(PetscInt));CHKERRQ(ierr); 637 } 638 PetscFunctionReturn(0); 639 } 640 641 642 #undef __FUNCT__ 643 #define __FUNCT__ "PetscProcessTree" 644 /*@ 645 PetscProcessTree - Prepares tree data to be displayed graphically 646 647 Not Collective 648 649 Input Parameters: 650 + n - number of values 651 . mask - indicates those entries in the tree, location 0 is always masked 652 - parentid - indicates the parent of each entry 653 654 Output Parameters: 655 + Nlevels - the number of levels 656 . Level - for each node tells its level 657 . Levelcnts - the number of nodes on each level 658 . Idbylevel - a list of ids on each of the levels, first level followed by second etc 659 - Column - for each id tells its column index 660 661 Level: intermediate 662 663 664 .seealso: PetscSortReal(), PetscSortIntWithPermutation() 665 @*/ 666 PetscErrorCode PetscProcessTree(PetscInt n,const PetscBool mask[],const PetscInt parentid[],PetscInt *Nlevels,PetscInt **Level,PetscInt **Levelcnt,PetscInt **Idbylevel,PetscInt **Column) 667 { 668 PetscInt i,j,cnt,nmask = 0,nlevels = 0,*level,*levelcnt,levelmax = 0,*workid,*workparentid,tcnt = 0,*idbylevel,*column; 669 PetscErrorCode ierr; 670 PetscBool done = PETSC_FALSE; 671 672 PetscFunctionBegin; 673 if (!mask[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mask of 0th location must be set"); 674 for (i=0; i<n; i++) { 675 if (mask[i]) continue; 676 if (parentid[i] == i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Node labeled as own parent"); 677 if (parentid[i] && mask[parentid[i]]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Parent is masked"); 678 } 679 680 for (i=0; i<n; i++) { 681 if (!mask[i]) nmask++; 682 } 683 684 /* determine the level in the tree of each node */ 685 ierr = PetscMalloc(n*sizeof(PetscInt),&level);CHKERRQ(ierr); 686 ierr = PetscMemzero(level,n*sizeof(PetscInt));CHKERRQ(ierr); 687 level[0] = 1; 688 while (!done) { 689 done = PETSC_TRUE; 690 for (i=0; i<n; i++) { 691 if (mask[i]) continue; 692 if (!level[i] && level[parentid[i]]) level[i] = level[parentid[i]] + 1; 693 else if (!level[i]) done = PETSC_FALSE; 694 } 695 } 696 for (i=0; i<n; i++) { 697 level[i]--; 698 nlevels = PetscMax(nlevels,level[i]); 699 } 700 701 /* count the number of nodes on each level and its max */ 702 ierr = PetscMalloc(nlevels*sizeof(PetscInt),&levelcnt);CHKERRQ(ierr); 703 ierr = PetscMemzero(levelcnt,nlevels*sizeof(PetscInt));CHKERRQ(ierr); 704 for (i=0; i<n; i++) { 705 if (mask[i]) continue; 706 levelcnt[level[i]-1]++; 707 } 708 for (i=0; i<nlevels;i++) { 709 levelmax = PetscMax(levelmax,levelcnt[i]); 710 } 711 712 /* for each level sort the ids by the parent id */ 713 ierr = PetscMalloc2(levelmax,PetscInt,&workid,levelmax,PetscInt,&workparentid);CHKERRQ(ierr); 714 ierr = PetscMalloc(nmask*sizeof(PetscInt),&idbylevel);CHKERRQ(ierr); 715 for (j=1; j<=nlevels;j++) { 716 cnt = 0; 717 for (i=0; i<n; i++) { 718 if (mask[i]) continue; 719 if (level[i] != j) continue; 720 workid[cnt] = i; 721 workparentid[cnt++] = parentid[i]; 722 } 723 /* PetscIntView(cnt,workparentid,0); 724 PetscIntView(cnt,workid,0); 725 ierr = PetscSortIntWithArray(cnt,workparentid,workid);CHKERRQ(ierr); 726 PetscIntView(cnt,workparentid,0); 727 PetscIntView(cnt,workid,0);*/ 728 ierr = PetscMemcpy(idbylevel+tcnt,workid,cnt*sizeof(PetscInt));CHKERRQ(ierr); 729 tcnt += cnt; 730 } 731 if (tcnt != nmask) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent count of unmasked nodes"); 732 ierr = PetscFree2(workid,workparentid);CHKERRQ(ierr); 733 734 /* for each node list its column */ 735 ierr = PetscMalloc(n*sizeof(PetscInt),&column);CHKERRQ(ierr); 736 cnt = 0; 737 for (j=0; j<nlevels; j++) { 738 for (i=0; i<levelcnt[j]; i++) { 739 column[idbylevel[cnt++]] = i; 740 } 741 } 742 743 *Nlevels = nlevels; 744 *Level = level; 745 *Levelcnt = levelcnt; 746 *Idbylevel = idbylevel; 747 *Column = column; 748 PetscFunctionReturn(0); 749 } 750