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 568 569 #undef __FUNCT__ 570 #define __FUNCT__ "PetscMergeIntArrayPair" 571 /*@ 572 PetscMergeIntArrayPair - Merges two SORTED integer arrays along with an additional array of integers. 573 The additional arrays are the same length as sorted arrays and are merged 574 in the order determined by the merging of the sorted pair. 575 576 Not Collective 577 578 Input Parameters: 579 + an - number of values in the first array 580 . aI - first sorted array of integers 581 . aJ - first additional array of integers 582 . bn - number of values in the second array 583 . bI - second array of integers 584 - bJ - second additional of integers 585 586 Output Parameters: 587 + n - number of values in the merged array (== an + bn) 588 . I - merged sorted array 589 - J - merged additional array 590 591 Level: intermediate 592 593 Concepts: merging^arrays 594 595 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray() 596 @*/ 597 PetscErrorCode PetscMergeIntArrayPair(PetscInt an,const PetscInt *aI, const PetscInt *aJ, PetscInt bn, const PetscInt *bI, const PetscInt *bJ, PetscInt *n, PetscInt **L, PetscInt **J) 598 { 599 PetscErrorCode ierr; 600 PetscInt n_, *L_ = *L, *J_= *J, ak, bk, k; 601 602 n_ = an + bn; 603 *n = n_; 604 if (!L_) { 605 ierr = PetscMalloc(n_*sizeof(PetscInt), L);CHKERRQ(ierr); 606 L_ = *L; 607 } 608 if (!J_) { 609 ierr = PetscMalloc(n_*sizeof(PetscInt), &J_);CHKERRQ(ierr); 610 J_ = *J; 611 } 612 k = ak = bk = 0; 613 while (ak < an && bk < bn) { 614 if (aI[ak] <= bI[bk]) { 615 L_[k] = aI[ak]; 616 J_[k] = aJ[ak]; 617 ++ak; 618 ++k; 619 } else { 620 L_[k] = bI[bk]; 621 J_[k] = bJ[bk]; 622 ++bk; 623 ++k; 624 } 625 } 626 if (ak < an) { 627 ierr = PetscMemcpy(L_+k,aI+ak,(an-ak)*sizeof(PetscInt));CHKERRQ(ierr); 628 ierr = PetscMemcpy(J_+k,aJ+ak,(an-ak)*sizeof(PetscInt));CHKERRQ(ierr); 629 k += (an-ak); 630 } 631 if (bk < bn) { 632 ierr = PetscMemcpy(L_+k,bI+bk,(bn-bk)*sizeof(PetscInt));CHKERRQ(ierr); 633 ierr = PetscMemcpy(J_+k,bJ+bk,(bn-bk)*sizeof(PetscInt));CHKERRQ(ierr); 634 } 635 PetscFunctionReturn(0); 636 } 637 638 639 #undef __FUNCT__ 640 #define __FUNCT__ "PetscProcessTree" 641 /*@ 642 PetscProcessTree - Prepares tree data to be displayed graphically 643 644 Not Collective 645 646 Input Parameters: 647 + n - number of values 648 . mask - indicates those entries in the tree, location 0 is always masked 649 - parentid - indicates the parent of each entry 650 651 Output Parameters: 652 + Nlevels - the number of levels 653 . Level - for each node tells its level 654 . Levelcnts - the number of nodes on each level 655 . Idbylevel - a list of ids on each of the levels, first level followed by second etc 656 - Column - for each id tells its column index 657 658 Level: intermediate 659 660 661 .seealso: PetscSortReal(), PetscSortIntWithPermutation() 662 @*/ 663 PetscErrorCode PetscProcessTree(PetscInt n,const PetscBool mask[],const PetscInt parentid[],PetscInt *Nlevels,PetscInt **Level,PetscInt **Levelcnt,PetscInt **Idbylevel,PetscInt **Column) 664 { 665 PetscInt i,j,cnt,nmask = 0,nlevels = 0,*level,*levelcnt,levelmax = 0,*workid,*workparentid,tcnt = 0,*idbylevel,*column; 666 PetscErrorCode ierr; 667 PetscBool done = PETSC_FALSE; 668 669 PetscFunctionBegin; 670 if (!mask[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mask of 0th location must be set"); 671 for (i=0; i<n; i++) { 672 if (mask[i]) continue; 673 if (parentid[i] == i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Node labeled as own parent"); 674 if (parentid[i] && mask[parentid[i]]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Parent is masked"); 675 } 676 677 for (i=0; i<n; i++) { 678 if (!mask[i]) nmask++; 679 } 680 681 /* determine the level in the tree of each node */ 682 ierr = PetscMalloc(n*sizeof(PetscInt),&level);CHKERRQ(ierr); 683 ierr = PetscMemzero(level,n*sizeof(PetscInt));CHKERRQ(ierr); 684 685 level[0] = 1; 686 while (!done) { 687 done = PETSC_TRUE; 688 for (i=0; i<n; i++) { 689 if (mask[i]) continue; 690 if (!level[i] && level[parentid[i]]) level[i] = level[parentid[i]] + 1; 691 else if (!level[i]) done = PETSC_FALSE; 692 } 693 } 694 for (i=0; i<n; i++) { 695 level[i]--; 696 nlevels = PetscMax(nlevels,level[i]); 697 } 698 699 /* count the number of nodes on each level and its max */ 700 ierr = PetscMalloc(nlevels*sizeof(PetscInt),&levelcnt);CHKERRQ(ierr); 701 ierr = PetscMemzero(levelcnt,nlevels*sizeof(PetscInt));CHKERRQ(ierr); 702 for (i=0; i<n; i++) { 703 if (mask[i]) continue; 704 levelcnt[level[i]-1]++; 705 } 706 for (i=0; i<nlevels;i++) levelmax = PetscMax(levelmax,levelcnt[i]); 707 708 /* for each level sort the ids by the parent id */ 709 ierr = PetscMalloc2(levelmax,PetscInt,&workid,levelmax,PetscInt,&workparentid);CHKERRQ(ierr); 710 ierr = PetscMalloc(nmask*sizeof(PetscInt),&idbylevel);CHKERRQ(ierr); 711 for (j=1; j<=nlevels;j++) { 712 cnt = 0; 713 for (i=0; i<n; i++) { 714 if (mask[i]) continue; 715 if (level[i] != j) continue; 716 workid[cnt] = i; 717 workparentid[cnt++] = parentid[i]; 718 } 719 /* PetscIntView(cnt,workparentid,0); 720 PetscIntView(cnt,workid,0); 721 ierr = PetscSortIntWithArray(cnt,workparentid,workid);CHKERRQ(ierr); 722 PetscIntView(cnt,workparentid,0); 723 PetscIntView(cnt,workid,0);*/ 724 ierr = PetscMemcpy(idbylevel+tcnt,workid,cnt*sizeof(PetscInt));CHKERRQ(ierr); 725 tcnt += cnt; 726 } 727 if (tcnt != nmask) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent count of unmasked nodes"); 728 ierr = PetscFree2(workid,workparentid);CHKERRQ(ierr); 729 730 /* for each node list its column */ 731 ierr = PetscMalloc(n*sizeof(PetscInt),&column);CHKERRQ(ierr); 732 cnt = 0; 733 for (j=0; j<nlevels; j++) { 734 for (i=0; i<levelcnt[j]; i++) { 735 column[idbylevel[cnt++]] = i; 736 } 737 } 738 739 *Nlevels = nlevels; 740 *Level = level; 741 *Levelcnt = levelcnt; 742 *Idbylevel = idbylevel; 743 *Column = column; 744 PetscFunctionReturn(0); 745 } 746