1 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*/ 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 PetscSortInt_Private(i,n-1); 86 PetscFunctionReturn(0); 87 } 88 89 #undef __FUNCT__ 90 #define __FUNCT__ "PetscSortRemoveDupsInt" 91 /*@ 92 PetscSortRemoveDupsInt - Sorts an array of integers in place in increasing order removes all duplicate entries 93 94 Not Collective 95 96 Input Parameters: 97 + n - number of values 98 - ii - array of integers 99 100 Output Parameter: 101 . n - number of non-redundant values 102 103 Level: intermediate 104 105 Concepts: sorting^ints 106 107 .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt() 108 @*/ 109 PetscErrorCode PetscSortRemoveDupsInt(PetscInt *n,PetscInt ii[]) 110 { 111 PetscErrorCode ierr; 112 PetscInt i,s = 0,N = *n, b = 0; 113 114 PetscFunctionBegin; 115 ierr = PetscSortInt(N,ii);CHKERRQ(ierr); 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 PetscFunctionReturn(0); 123 } 124 125 #undef __FUNCT__ 126 #define __FUNCT__ "PetscFindInt" 127 /*@ 128 PetscFindInt - Finds integer in a sorted array of integers 129 130 Not Collective 131 132 Input Parameters: 133 + key - the integer to locate 134 . n - number of values in the array 135 - ii - array of integers 136 137 Output Parameter: 138 . loc - the location if found, otherwise -(slot+1) where slot is the place the value would go 139 140 Level: intermediate 141 142 Concepts: sorting^ints 143 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; 149 150 PetscFunctionBegin; 151 PetscValidPointer(loc,4); 152 if (!n) {*loc = -1; PetscFunctionReturn(0);} 153 PetscValidPointer(ii,3); 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 PetscSortMPIInt_Private(i,n-1); 381 PetscFunctionReturn(0); 382 } 383 384 #undef __FUNCT__ 385 #define __FUNCT__ "PetscSortRemoveDupsMPIInt" 386 /*@ 387 PetscSortRemoveDupsMPIInt - Sorts an array of MPI integers in place in increasing order removes all duplicate entries 388 389 Not Collective 390 391 Input Parameters: 392 + n - number of values 393 - ii - array of integers 394 395 Output Parameter: 396 . n - number of non-redundant values 397 398 Level: intermediate 399 400 Concepts: sorting^ints 401 402 .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt() 403 @*/ 404 PetscErrorCode PetscSortRemoveDupsMPIInt(PetscInt *n,PetscMPIInt ii[]) 405 { 406 PetscErrorCode ierr; 407 PetscInt i,s = 0,N = *n, b = 0; 408 409 PetscFunctionBegin; 410 ierr = PetscSortMPIInt(N,ii);CHKERRQ(ierr); 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 PetscFunctionReturn(0); 418 } 419 420 #undef __FUNCT__ 421 #define __FUNCT__ "PetscSortMPIIntWithArray_Private" 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 { 429 PetscErrorCode ierr; 430 PetscMPIInt i,vl,last,tmp; 431 432 PetscFunctionBegin; 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 PetscFunctionReturn(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 ierr = PetscSortMPIIntWithArray_Private(v,V,last-1);CHKERRQ(ierr); 447 ierr = PetscSortMPIIntWithArray_Private(v+last+1,V+last+1,right-(last+1));CHKERRQ(ierr); 448 PetscFunctionReturn(0); 449 } 450 451 #undef __FUNCT__ 452 #define __FUNCT__ "PetscSortMPIIntWithArray" 453 /*@ 454 PetscSortMPIIntWithArray - Sorts an array of integers in place in increasing order; 455 changes a second array to match the sorted first array. 456 457 Not Collective 458 459 Input Parameters: 460 + n - number of values 461 . i - array of integers 462 - I - second array of integers 463 464 Level: intermediate 465 466 Concepts: sorting^ints with array 467 468 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt() 469 @*/ 470 PetscErrorCode PetscSortMPIIntWithArray(PetscMPIInt n,PetscMPIInt i[],PetscMPIInt Ii[]) 471 { 472 PetscErrorCode ierr; 473 PetscMPIInt j,k,tmp,ik; 474 475 PetscFunctionBegin; 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 ierr = PetscSortMPIIntWithArray_Private(i,Ii,n-1);CHKERRQ(ierr); 488 } 489 PetscFunctionReturn(0); 490 } 491 492 /* -----------------------------------------------------------------------*/ 493 #define SWAP2IntScalar(a,b,c,d,t,ts) {t=a;a=b;b=t;ts=c;c=d;d=ts;} 494 495 #undef __FUNCT__ 496 #define __FUNCT__ "PetscSortIntWithScalarArray_Private" 497 /* 498 Modified from PetscSortIntWithArray_Private(). 499 */ 500 static PetscErrorCode PetscSortIntWithScalarArray_Private(PetscInt *v,PetscScalar *V,PetscInt right) 501 { 502 PetscErrorCode ierr; 503 PetscInt i,vl,last,tmp; 504 PetscScalar stmp; 505 506 PetscFunctionBegin; 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 PetscFunctionReturn(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 ierr = PetscSortIntWithScalarArray_Private(v,V,last-1);CHKERRQ(ierr); 521 ierr = PetscSortIntWithScalarArray_Private(v+last+1,V+last+1,right-(last+1));CHKERRQ(ierr); 522 PetscFunctionReturn(0); 523 } 524 525 #undef __FUNCT__ 526 #define __FUNCT__ "PetscSortIntWithScalarArray" 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. 530 531 Not Collective 532 533 Input Parameters: 534 + n - number of values 535 . i - array of integers 536 - I - second array of scalars 537 538 Level: intermediate 539 540 Concepts: sorting^ints with array 541 542 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray() 543 @*/ 544 PetscErrorCode PetscSortIntWithScalarArray(PetscInt n,PetscInt i[],PetscScalar Ii[]) 545 { 546 PetscErrorCode ierr; 547 PetscInt j,k,tmp,ik; 548 PetscScalar stmp; 549 550 PetscFunctionBegin; 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 ierr = PetscSortIntWithScalarArray_Private(i,Ii,n-1);CHKERRQ(ierr); 563 } 564 PetscFunctionReturn(0); 565 } 566 567 #undef __FUNCT__ 568 #define __FUNCT__ "PetscMergeIntArrayPair" 569 /*@ 570 PetscMergeIntArrayPair - Merges two SORTED integer arrays along with an additional array of integers. 571 The additional arrays are the same length as sorted arrays and are merged 572 in the order determined by the merging of the sorted pair. 573 574 Not Collective 575 576 Input Parameters: 577 + an - number of values in the first array 578 . aI - first sorted array of integers 579 . aJ - first additional array of integers 580 . bn - number of values in the second array 581 . bI - second array of integers 582 - bJ - second additional of integers 583 584 Output Parameters: 585 + n - number of values in the merged array (== an + bn) 586 . L - merged sorted array 587 - J - merged additional array 588 589 Level: intermediate 590 591 Concepts: merging^arrays 592 593 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray() 594 @*/ 595 PetscErrorCode PetscMergeIntArrayPair(PetscInt an,const PetscInt *aI, const PetscInt *aJ, PetscInt bn, const PetscInt *bI, const PetscInt *bJ, PetscInt *n, PetscInt **L, PetscInt **J) 596 { 597 PetscErrorCode ierr; 598 PetscInt n_, *L_ = *L, *J_= *J, ak, bk, k; 599 600 PetscFunctionBegin; 601 n_ = an + bn; 602 *n = n_; 603 if (!L_) { 604 ierr = PetscMalloc1(n_, L);CHKERRQ(ierr); 605 L_ = *L; 606 } 607 if (!J_) { 608 ierr = PetscMalloc1(n_, &J_);CHKERRQ(ierr); 609 J_ = *J; 610 } 611 k = ak = bk = 0; 612 while (ak < an && bk < bn) { 613 if (aI[ak] <= bI[bk]) { 614 L_[k] = aI[ak]; 615 J_[k] = aJ[ak]; 616 ++ak; 617 ++k; 618 } else { 619 L_[k] = bI[bk]; 620 J_[k] = bJ[bk]; 621 ++bk; 622 ++k; 623 } 624 } 625 if (ak < an) { 626 ierr = PetscMemcpy(L_+k,aI+ak,(an-ak)*sizeof(PetscInt));CHKERRQ(ierr); 627 ierr = PetscMemcpy(J_+k,aJ+ak,(an-ak)*sizeof(PetscInt));CHKERRQ(ierr); 628 k += (an-ak); 629 } 630 if (bk < bn) { 631 ierr = PetscMemcpy(L_+k,bI+bk,(bn-bk)*sizeof(PetscInt));CHKERRQ(ierr); 632 ierr = PetscMemcpy(J_+k,bJ+bk,(bn-bk)*sizeof(PetscInt));CHKERRQ(ierr); 633 } 634 PetscFunctionReturn(0); 635 } 636 637 638 #undef __FUNCT__ 639 #define __FUNCT__ "PetscProcessTree" 640 /*@ 641 PetscProcessTree - Prepares tree data to be displayed graphically 642 643 Not Collective 644 645 Input Parameters: 646 + n - number of values 647 . mask - indicates those entries in the tree, location 0 is always masked 648 - parentid - indicates the parent of each entry 649 650 Output Parameters: 651 + Nlevels - the number of levels 652 . Level - for each node tells its level 653 . Levelcnts - the number of nodes on each level 654 . Idbylevel - a list of ids on each of the levels, first level followed by second etc 655 - Column - for each id tells its column index 656 657 Level: intermediate 658 659 660 .seealso: PetscSortReal(), PetscSortIntWithPermutation() 661 @*/ 662 PetscErrorCode PetscProcessTree(PetscInt n,const PetscBool mask[],const PetscInt parentid[],PetscInt *Nlevels,PetscInt **Level,PetscInt **Levelcnt,PetscInt **Idbylevel,PetscInt **Column) 663 { 664 PetscInt i,j,cnt,nmask = 0,nlevels = 0,*level,*levelcnt,levelmax = 0,*workid,*workparentid,tcnt = 0,*idbylevel,*column; 665 PetscErrorCode ierr; 666 PetscBool done = PETSC_FALSE; 667 668 PetscFunctionBegin; 669 if (!mask[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mask of 0th location must be set"); 670 for (i=0; i<n; i++) { 671 if (mask[i]) continue; 672 if (parentid[i] == i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Node labeled as own parent"); 673 if (parentid[i] && mask[parentid[i]]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Parent is masked"); 674 } 675 676 for (i=0; i<n; i++) { 677 if (!mask[i]) nmask++; 678 } 679 680 /* determine the level in the tree of each node */ 681 ierr = PetscCalloc1(n,&level);CHKERRQ(ierr); 682 683 level[0] = 1; 684 while (!done) { 685 done = PETSC_TRUE; 686 for (i=0; i<n; i++) { 687 if (mask[i]) continue; 688 if (!level[i] && level[parentid[i]]) level[i] = level[parentid[i]] + 1; 689 else if (!level[i]) done = PETSC_FALSE; 690 } 691 } 692 for (i=0; i<n; i++) { 693 level[i]--; 694 nlevels = PetscMax(nlevels,level[i]); 695 } 696 697 /* count the number of nodes on each level and its max */ 698 ierr = PetscCalloc1(nlevels,&levelcnt);CHKERRQ(ierr); 699 for (i=0; i<n; i++) { 700 if (mask[i]) continue; 701 levelcnt[level[i]-1]++; 702 } 703 for (i=0; i<nlevels;i++) levelmax = PetscMax(levelmax,levelcnt[i]); 704 705 /* for each level sort the ids by the parent id */ 706 ierr = PetscMalloc2(levelmax,&workid,levelmax,&workparentid);CHKERRQ(ierr); 707 ierr = PetscMalloc1(nmask,&idbylevel);CHKERRQ(ierr); 708 for (j=1; j<=nlevels;j++) { 709 cnt = 0; 710 for (i=0; i<n; i++) { 711 if (mask[i]) continue; 712 if (level[i] != j) continue; 713 workid[cnt] = i; 714 workparentid[cnt++] = parentid[i]; 715 } 716 /* PetscIntView(cnt,workparentid,0); 717 PetscIntView(cnt,workid,0); 718 ierr = PetscSortIntWithArray(cnt,workparentid,workid);CHKERRQ(ierr); 719 PetscIntView(cnt,workparentid,0); 720 PetscIntView(cnt,workid,0);*/ 721 ierr = PetscMemcpy(idbylevel+tcnt,workid,cnt*sizeof(PetscInt));CHKERRQ(ierr); 722 tcnt += cnt; 723 } 724 if (tcnt != nmask) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent count of unmasked nodes"); 725 ierr = PetscFree2(workid,workparentid);CHKERRQ(ierr); 726 727 /* for each node list its column */ 728 ierr = PetscMalloc1(n,&column);CHKERRQ(ierr); 729 cnt = 0; 730 for (j=0; j<nlevels; j++) { 731 for (i=0; i<levelcnt[j]; i++) { 732 column[idbylevel[cnt++]] = i; 733 } 734 } 735 736 *Nlevels = nlevels; 737 *Level = level; 738 *Levelcnt = levelcnt; 739 *Idbylevel = idbylevel; 740 *Column = column; 741 PetscFunctionReturn(0); 742 } 743