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 #undef __FUNCT__ 569 #define __FUNCT__ "PetscMergeIntArray" 570 /*@ 571 PetscMergeIntArray - Merges two SORTED integer arrays, removes duplicate elements. 572 573 Not Collective 574 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 580 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 584 585 Level: intermediate 586 587 Concepts: merging^arrays 588 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 { 593 PetscErrorCode ierr; 594 PetscInt *L_ = *L, ak, bk, k; 595 596 if (!L_) { 597 ierr = PetscMalloc1(an+bn, L);CHKERRQ(ierr); 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 ierr = PetscMemcpy(L_+k,aI+ak,(an-ak)*sizeof(PetscInt));CHKERRQ(ierr); 619 k += (an-ak); 620 } 621 if (bk < bn) { 622 ierr = PetscMemcpy(L_+k,bI+bk,(bn-bk)*sizeof(PetscInt));CHKERRQ(ierr); 623 k += (bn-bk); 624 } 625 *n = k; 626 PetscFunctionReturn(0); 627 } 628 629 #undef __FUNCT__ 630 #define __FUNCT__ "PetscMergeIntArrayPair" 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. 635 636 Not Collective 637 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 645 646 Output Parameters: 647 + n - number of values in the merged array (== an + bn) 648 . I - merged sorted array 649 - J - merged additional array 650 651 Level: intermediate 652 653 Concepts: merging^arrays 654 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 { 659 PetscErrorCode ierr; 660 PetscInt n_, *L_ = *L, *J_= *J, ak, bk, k; 661 662 n_ = an + bn; 663 *n = n_; 664 if (!L_) { 665 ierr = PetscMalloc1(n_, L);CHKERRQ(ierr); 666 L_ = *L; 667 } 668 if (!J_) { 669 ierr = PetscMalloc1(n_, &J_);CHKERRQ(ierr); 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 ierr = PetscMemcpy(L_+k,aI+ak,(an-ak)*sizeof(PetscInt));CHKERRQ(ierr); 688 ierr = PetscMemcpy(J_+k,aJ+ak,(an-ak)*sizeof(PetscInt));CHKERRQ(ierr); 689 k += (an-ak); 690 } 691 if (bk < bn) { 692 ierr = PetscMemcpy(L_+k,bI+bk,(bn-bk)*sizeof(PetscInt));CHKERRQ(ierr); 693 ierr = PetscMemcpy(J_+k,bJ+bk,(bn-bk)*sizeof(PetscInt));CHKERRQ(ierr); 694 } 695 PetscFunctionReturn(0); 696 } 697 698 699 #undef __FUNCT__ 700 #define __FUNCT__ "PetscProcessTree" 701 /*@ 702 PetscProcessTree - Prepares tree data to be displayed graphically 703 704 Not Collective 705 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 710 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 717 718 Level: intermediate 719 720 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; 726 PetscErrorCode ierr; 727 PetscBool done = PETSC_FALSE; 728 729 PetscFunctionBegin; 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 } 736 737 for (i=0; i<n; i++) { 738 if (!mask[i]) nmask++; 739 } 740 741 /* determine the level in the tree of each node */ 742 ierr = PetscCalloc1(n,&level);CHKERRQ(ierr); 743 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 } 757 758 /* count the number of nodes on each level and its max */ 759 ierr = PetscCalloc1(nlevels,&levelcnt);CHKERRQ(ierr); 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]); 765 766 /* for each level sort the ids by the parent id */ 767 ierr = PetscMalloc2(levelmax,&workid,levelmax,&workparentid);CHKERRQ(ierr); 768 ierr = PetscMalloc1(nmask,&idbylevel);CHKERRQ(ierr); 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 ierr = PetscSortIntWithArray(cnt,workparentid,workid);CHKERRQ(ierr); 780 PetscIntView(cnt,workparentid,0); 781 PetscIntView(cnt,workid,0);*/ 782 ierr = PetscMemcpy(idbylevel+tcnt,workid,cnt*sizeof(PetscInt));CHKERRQ(ierr); 783 tcnt += cnt; 784 } 785 if (tcnt != nmask) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent count of unmasked nodes"); 786 ierr = PetscFree2(workid,workparentid);CHKERRQ(ierr); 787 788 /* for each node list its column */ 789 ierr = PetscMalloc1(n,&column);CHKERRQ(ierr); 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 } 796 797 *Nlevels = nlevels; 798 *Level = level; 799 *Levelcnt = levelcnt; 800 *Idbylevel = idbylevel; 801 *Column = column; 802 PetscFunctionReturn(0); 803 } 804