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(ii,3); 153 PetscValidPointer(loc,4); 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 PetscFunctionReturn(0); 161 } 162 163 164 /* -----------------------------------------------------------------------*/ 165 #define SWAP2(a,b,c,d,t) {t=a;a=b;b=t;t=c;c=d;d=t;} 166 167 #undef __FUNCT__ 168 #define __FUNCT__ "PetscSortIntWithArray_Private" 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 { 176 PetscErrorCode ierr; 177 PetscInt i,vl,last,tmp; 178 179 PetscFunctionBegin; 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 PetscFunctionReturn(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 ierr = PetscSortIntWithArray_Private(v,V,last-1);CHKERRQ(ierr); 194 ierr = PetscSortIntWithArray_Private(v+last+1,V+last+1,right-(last+1));CHKERRQ(ierr); 195 PetscFunctionReturn(0); 196 } 197 198 #undef __FUNCT__ 199 #define __FUNCT__ "PetscSortIntWithArray" 200 /*@ 201 PetscSortIntWithArray - Sorts an array of integers in place in increasing order; 202 changes a second array to match the sorted first array. 203 204 Not Collective 205 206 Input Parameters: 207 + n - number of values 208 . i - array of integers 209 - I - second array of integers 210 211 Level: intermediate 212 213 Concepts: sorting^ints with array 214 215 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt() 216 @*/ 217 PetscErrorCode PetscSortIntWithArray(PetscInt n,PetscInt i[],PetscInt Ii[]) 218 { 219 PetscErrorCode ierr; 220 PetscInt j,k,tmp,ik; 221 222 PetscFunctionBegin; 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 ierr = PetscSortIntWithArray_Private(i,Ii,n-1);CHKERRQ(ierr); 235 } 236 PetscFunctionReturn(0); 237 } 238 239 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;} 241 242 #undef __FUNCT__ 243 #define __FUNCT__ "PetscSortIntWithArrayPair_Private" 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 { 251 PetscErrorCode ierr; 252 PetscInt i,vl,last,tmp; 253 254 PetscFunctionBegin; 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 PetscFunctionReturn(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 ierr = PetscSortIntWithArrayPair_Private(L,J,K,last-1);CHKERRQ(ierr); 269 ierr = PetscSortIntWithArrayPair_Private(L+last+1,J+last+1,K+last+1,right-(last+1));CHKERRQ(ierr); 270 PetscFunctionReturn(0); 271 } 272 273 #undef __FUNCT__ 274 #define __FUNCT__ "PetscSortIntWithArrayPair" 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. 278 279 Not Collective 280 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) 286 287 Level: intermediate 288 289 Concepts: sorting^ints with array pair 290 291 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortIntWithArray() 292 @*/ 293 PetscErrorCode PetscSortIntWithArrayPair(PetscInt n,PetscInt *L,PetscInt *J, PetscInt *K) 294 { 295 PetscErrorCode ierr; 296 PetscInt j,k,tmp,ik; 297 298 PetscFunctionBegin; 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 ierr = PetscSortIntWithArrayPair_Private(L,J,K,n-1);CHKERRQ(ierr); 311 } 312 PetscFunctionReturn(0); 313 } 314 315 #undef __FUNCT__ 316 #define __FUNCT__ "PetscSortMPIInt_Private" 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; 326 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 } 346 347 #undef __FUNCT__ 348 #define __FUNCT__ "PetscSortMPIInt" 349 /*@ 350 PetscSortMPIInt - Sorts an array of MPI integers in place in increasing order. 351 352 Not Collective 353 354 Input Parameters: 355 + n - number of values 356 - i - array of integers 357 358 Level: intermediate 359 360 Concepts: sorting^ints 361 362 .seealso: PetscSortReal(), PetscSortIntWithPermutation() 363 @*/ 364 PetscErrorCode PetscSortMPIInt(PetscInt n,PetscMPIInt i[]) 365 { 366 PetscInt j,k; 367 PetscMPIInt tmp,ik; 368 369 PetscFunctionBegin; 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 { 381 PetscSortMPIInt_Private(i,n-1); 382 } 383 PetscFunctionReturn(0); 384 } 385 386 #undef __FUNCT__ 387 #define __FUNCT__ "PetscSortRemoveDupsMPIInt" 388 /*@ 389 PetscSortRemoveDupsMPIInt - Sorts an array of MPI integers in place in increasing order removes all duplicate entries 390 391 Not Collective 392 393 Input Parameters: 394 + n - number of values 395 - ii - array of integers 396 397 Output Parameter: 398 . n - number of non-redundant values 399 400 Level: intermediate 401 402 Concepts: sorting^ints 403 404 .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt() 405 @*/ 406 PetscErrorCode PetscSortRemoveDupsMPIInt(PetscInt *n,PetscMPIInt ii[]) 407 { 408 PetscErrorCode ierr; 409 PetscInt i,s = 0,N = *n, b = 0; 410 411 PetscFunctionBegin; 412 ierr = PetscSortMPIInt(N,ii);CHKERRQ(ierr); 413 for (i=0; i<N-1; i++) { 414 if (ii[b+s+1] != ii[b]) {ii[b+1] = ii[b+s+1]; b++;} 415 else s++; 416 } 417 *n = N - s; 418 PetscFunctionReturn(0); 419 } 420 421 #undef __FUNCT__ 422 #define __FUNCT__ "PetscSortMPIIntWithArray_Private" 423 /* 424 A simple version of quicksort; taken from Kernighan and Ritchie, page 87. 425 Assumes 0 origin for v, number of elements = right+1 (right is index of 426 right-most member). 427 */ 428 static PetscErrorCode PetscSortMPIIntWithArray_Private(PetscMPIInt *v,PetscMPIInt *V,PetscMPIInt right) 429 { 430 PetscErrorCode ierr; 431 PetscMPIInt i,vl,last,tmp; 432 433 PetscFunctionBegin; 434 if (right <= 1) { 435 if (right == 1) { 436 if (v[0] > v[1]) SWAP2(v[0],v[1],V[0],V[1],tmp); 437 } 438 PetscFunctionReturn(0); 439 } 440 SWAP2(v[0],v[right/2],V[0],V[right/2],tmp); 441 vl = v[0]; 442 last = 0; 443 for (i=1; i<=right; i++) { 444 if (v[i] < vl) {last++; SWAP2(v[last],v[i],V[last],V[i],tmp);} 445 } 446 SWAP2(v[0],v[last],V[0],V[last],tmp); 447 ierr = PetscSortMPIIntWithArray_Private(v,V,last-1);CHKERRQ(ierr); 448 ierr = PetscSortMPIIntWithArray_Private(v+last+1,V+last+1,right-(last+1));CHKERRQ(ierr); 449 PetscFunctionReturn(0); 450 } 451 452 #undef __FUNCT__ 453 #define __FUNCT__ "PetscSortMPIIntWithArray" 454 /*@ 455 PetscSortMPIIntWithArray - Sorts an array of integers in place in increasing order; 456 changes a second array to match the sorted first array. 457 458 Not Collective 459 460 Input Parameters: 461 + n - number of values 462 . i - array of integers 463 - I - second array of integers 464 465 Level: intermediate 466 467 Concepts: sorting^ints with array 468 469 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt() 470 @*/ 471 PetscErrorCode PetscSortMPIIntWithArray(PetscMPIInt n,PetscMPIInt i[],PetscMPIInt Ii[]) 472 { 473 PetscErrorCode ierr; 474 PetscMPIInt j,k,tmp,ik; 475 476 PetscFunctionBegin; 477 if (n<8) { 478 for (k=0; k<n; k++) { 479 ik = i[k]; 480 for (j=k+1; j<n; j++) { 481 if (ik > i[j]) { 482 SWAP2(i[k],i[j],Ii[k],Ii[j],tmp); 483 ik = i[k]; 484 } 485 } 486 } 487 } else { 488 ierr = PetscSortMPIIntWithArray_Private(i,Ii,n-1);CHKERRQ(ierr); 489 } 490 PetscFunctionReturn(0); 491 } 492 493 /* -----------------------------------------------------------------------*/ 494 #define SWAP2IntScalar(a,b,c,d,t,ts) {t=a;a=b;b=t;ts=c;c=d;d=ts;} 495 496 #undef __FUNCT__ 497 #define __FUNCT__ "PetscSortIntWithScalarArray_Private" 498 /* 499 Modified from PetscSortIntWithArray_Private(). 500 */ 501 static PetscErrorCode PetscSortIntWithScalarArray_Private(PetscInt *v,PetscScalar *V,PetscInt right) 502 { 503 PetscErrorCode ierr; 504 PetscInt i,vl,last,tmp; 505 PetscScalar stmp; 506 507 PetscFunctionBegin; 508 if (right <= 1) { 509 if (right == 1) { 510 if (v[0] > v[1]) SWAP2IntScalar(v[0],v[1],V[0],V[1],tmp,stmp); 511 } 512 PetscFunctionReturn(0); 513 } 514 SWAP2IntScalar(v[0],v[right/2],V[0],V[right/2],tmp,stmp); 515 vl = v[0]; 516 last = 0; 517 for (i=1; i<=right; i++) { 518 if (v[i] < vl) {last++; SWAP2IntScalar(v[last],v[i],V[last],V[i],tmp,stmp);} 519 } 520 SWAP2IntScalar(v[0],v[last],V[0],V[last],tmp,stmp); 521 ierr = PetscSortIntWithScalarArray_Private(v,V,last-1);CHKERRQ(ierr); 522 ierr = PetscSortIntWithScalarArray_Private(v+last+1,V+last+1,right-(last+1));CHKERRQ(ierr); 523 PetscFunctionReturn(0); 524 } 525 526 #undef __FUNCT__ 527 #define __FUNCT__ "PetscSortIntWithScalarArray" 528 /*@ 529 PetscSortIntWithScalarArray - Sorts an array of integers in place in increasing order; 530 changes a second SCALAR array to match the sorted first INTEGER array. 531 532 Not Collective 533 534 Input Parameters: 535 + n - number of values 536 . i - array of integers 537 - I - second array of scalars 538 539 Level: intermediate 540 541 Concepts: sorting^ints with array 542 543 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray() 544 @*/ 545 PetscErrorCode PetscSortIntWithScalarArray(PetscInt n,PetscInt i[],PetscScalar Ii[]) 546 { 547 PetscErrorCode ierr; 548 PetscInt j,k,tmp,ik; 549 PetscScalar stmp; 550 551 PetscFunctionBegin; 552 if (n<8) { 553 for (k=0; k<n; k++) { 554 ik = i[k]; 555 for (j=k+1; j<n; j++) { 556 if (ik > i[j]) { 557 SWAP2IntScalar(i[k],i[j],Ii[k],Ii[j],tmp,stmp); 558 ik = i[k]; 559 } 560 } 561 } 562 } else { 563 ierr = PetscSortIntWithScalarArray_Private(i,Ii,n-1);CHKERRQ(ierr); 564 } 565 PetscFunctionReturn(0); 566 } 567 568 569 570 #undef __FUNCT__ 571 #define __FUNCT__ "PetscMergeIntArrayPair" 572 /*@ 573 PetscMergeIntArrayPair - Merges two SORTED integer arrays along with an additional array of integers. 574 The additional arrays are the same length as sorted arrays and are merged 575 in the order determined by the merging of the sorted pair. 576 577 Not Collective 578 579 Input Parameters: 580 + an - number of values in the first array 581 . aI - first sorted array of integers 582 . aJ - first additional array of integers 583 . bn - number of values in the second array 584 . bI - second array of integers 585 - bJ - second additional of integers 586 587 Output Parameters: 588 + n - number of values in the merged array (== an + bn) 589 . I - merged sorted array 590 - J - merged additional array 591 592 Level: intermediate 593 594 Concepts: merging^arrays 595 596 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray() 597 @*/ 598 PetscErrorCode PetscMergeIntArrayPair(PetscInt an,const PetscInt *aI, const PetscInt *aJ, PetscInt bn, const PetscInt *bI, const PetscInt *bJ, PetscInt *n, PetscInt **L, PetscInt **J) 599 { 600 PetscErrorCode ierr; 601 PetscInt n_, *L_ = *L, *J_= *J, ak, bk, k; 602 603 n_ = an + bn; 604 *n = n_; 605 if (!L_) { 606 ierr = PetscMalloc(n_*sizeof(PetscInt), L); CHKERRQ(ierr); 607 L_ = *L; 608 } 609 if (!J_){ 610 ierr = PetscMalloc(n_*sizeof(PetscInt), &J_); CHKERRQ(ierr); 611 J_ = *J; 612 } 613 k = ak = bk = 0; 614 while(ak < an && bk < bn) { 615 if (aI[ak] <= bI[bk]) { 616 L_[k] = aI[ak]; 617 J_[k] = aJ[ak]; 618 ++ak; 619 ++k; 620 } 621 else { 622 L_[k] = bI[bk]; 623 J_[k] = bJ[bk]; 624 ++bk; 625 ++k; 626 } 627 } 628 if (ak < an) { 629 ierr = PetscMemcpy(L_+k,aI+ak,(an-ak)*sizeof(PetscInt)); CHKERRQ(ierr); 630 ierr = PetscMemcpy(J_+k,aJ+ak,(an-ak)*sizeof(PetscInt)); CHKERRQ(ierr); 631 k += (an-ak); 632 } 633 if (bk < bn) { 634 ierr = PetscMemcpy(L_+k,bI+bk,(bn-bk)*sizeof(PetscInt)); CHKERRQ(ierr); 635 ierr = PetscMemcpy(J_+k,bJ+bk,(bn-bk)*sizeof(PetscInt)); CHKERRQ(ierr); 636 k += (bn-bk); 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 681 for (i=0; i<n; i++) { 682 if (!mask[i]) nmask++; 683 } 684 685 /* determine the level in the tree of each node */ 686 ierr = PetscMalloc(n*sizeof(PetscInt),&level);CHKERRQ(ierr); 687 ierr = PetscMemzero(level,n*sizeof(PetscInt));CHKERRQ(ierr); 688 level[0] = 1; 689 while (!done) { 690 done = PETSC_TRUE; 691 for (i=0; i<n; i++) { 692 if (mask[i]) continue; 693 if (!level[i] && level[parentid[i]]) level[i] = level[parentid[i]] + 1; 694 else if (!level[i]) done = PETSC_FALSE; 695 } 696 } 697 for (i=0; i<n; i++) { 698 level[i]--; 699 nlevels = PetscMax(nlevels,level[i]); 700 } 701 702 /* count the number of nodes on each level and its max */ 703 ierr = PetscMalloc(nlevels*sizeof(PetscInt),&levelcnt);CHKERRQ(ierr); 704 ierr = PetscMemzero(levelcnt,nlevels*sizeof(PetscInt));CHKERRQ(ierr); 705 for (i=0; i<n; i++) { 706 if (mask[i]) continue; 707 levelcnt[level[i]-1]++; 708 } 709 for (i=0; i<nlevels;i++) { 710 levelmax = PetscMax(levelmax,levelcnt[i]); 711 } 712 713 /* for each level sort the ids by the parent id */ 714 ierr = PetscMalloc2(levelmax,PetscInt,&workid,levelmax,PetscInt,&workparentid);CHKERRQ(ierr); 715 ierr = PetscMalloc(nmask*sizeof(PetscInt),&idbylevel);CHKERRQ(ierr); 716 for (j=1; j<=nlevels;j++) { 717 cnt = 0; 718 for (i=0; i<n; i++) { 719 if (mask[i]) continue; 720 if (level[i] != j) continue; 721 workid[cnt] = i; 722 workparentid[cnt++] = parentid[i]; 723 } 724 /* PetscIntView(cnt,workparentid,0); 725 PetscIntView(cnt,workid,0); 726 ierr = PetscSortIntWithArray(cnt,workparentid,workid);CHKERRQ(ierr); 727 PetscIntView(cnt,workparentid,0); 728 PetscIntView(cnt,workid,0);*/ 729 ierr = PetscMemcpy(idbylevel+tcnt,workid,cnt*sizeof(PetscInt));CHKERRQ(ierr); 730 tcnt += cnt; 731 } 732 if (tcnt != nmask) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent count of unmasked nodes"); 733 ierr = PetscFree2(workid,workparentid);CHKERRQ(ierr); 734 735 /* for each node list its column */ 736 ierr = PetscMalloc(n*sizeof(PetscInt),&column);CHKERRQ(ierr); 737 cnt = 0; 738 for (j=0; j<nlevels; j++) { 739 for (i=0; i<levelcnt[j]; i++) { 740 column[idbylevel[cnt++]] = i; 741 } 742 } 743 744 *Nlevels = nlevels; 745 *Level = level; 746 *Levelcnt = levelcnt; 747 *Idbylevel = idbylevel; 748 *Column = column; 749 PetscFunctionReturn(0); 750 } 751