1 2 /* 3 This file contains routines for sorting integers. Values are sorted in place. 4 One can use src/sys/examples/tests/ex52.c for benchmarking. 5 */ 6 #include <petsc/private/petscimpl.h> /*I "petscsys.h" I*/ 7 #include <petsc/private/hashseti.h> 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 /* Swap one, two or three pairs. Each pair can have its own type */ 21 #define SWAP1(a,b,t1) do {t1=a;a=b;b=t1;} while(0) 22 #define SWAP2(a,b,c,d,t1,t2) do {t1=a;a=b;b=t1; t2=c;c=d;d=t2;} while(0) 23 #define SWAP3(a,b,c,d,e,f,t1,t2,t3) do {t1=a;a=b;b=t1; t2=c;c=d;d=t2; t3=e;e=f;f=t3;} while(0) 24 25 /* Swap a & b, *c & *d. c, d, t2 are pointers to a type of size <siz> */ 26 #define SWAP2Data(a,b,c,d,t1,t2,siz) \ 27 do { \ 28 PetscErrorCode ierr; \ 29 t1=a; a=b; b=t1; \ 30 ierr = PetscMemcpy(t2,c,siz);CHKERRQ(ierr); \ 31 ierr = PetscMemcpy(c,d,siz);CHKERRQ(ierr); \ 32 ierr = PetscMemcpy(d,t2,siz);CHKERRQ(ierr); \ 33 } while(0) 34 35 /* 36 Partition X[lo,hi] into two parts: X[lo,l) <= pivot; X[r,hi] > pivot 37 38 Input Parameters: 39 + X - array to partition 40 . pivot - a pivot of X[] 41 . t1 - temp variable for X 42 - lo,hi - lower and upper bound of the array 43 44 Output Parameters: 45 + l,r - of type PetscInt 46 47 Notes: 48 The TwoWayPartition2/3 variants also partition other arrays along with X. 49 These arrays can have different types, so they provide their own temp t2,t3 50 */ 51 #define TwoWayPartition1(X,pivot,t1,lo,hi,l,r) \ 52 do { \ 53 l = lo; \ 54 r = hi; \ 55 while(1) { \ 56 while (X[l] < pivot) l++; \ 57 while (X[r] > pivot) r--; \ 58 if (l >= r) {r++; break;} \ 59 SWAP1(X[l],X[r],t1); \ 60 l++; \ 61 r--; \ 62 } \ 63 } while(0) 64 65 /* 66 Partition X[lo,hi] into two parts: X[lo,l) >= pivot; X[r,hi] < pivot 67 68 Input Parameters: 69 + X - array to partition 70 . pivot - a pivot of X[] 71 . t1 - temp variable for X 72 - lo,hi - lower and upper bound of the array 73 74 Output Parameters: 75 + l,r - of type PetscInt 76 77 Notes: 78 The TwoWayPartition2/3 variants also partition other arrays along with X. 79 These arrays can have different types, so they provide their own temp t2,t3 80 */ 81 #define TwoWayPartitionReverse1(X,pivot,t1,lo,hi,l,r) \ 82 do { \ 83 l = lo; \ 84 r = hi; \ 85 while(1) { \ 86 while (X[l] > pivot) l++; \ 87 while (X[r] < pivot) r--; \ 88 if (l >= r) {r++; break;} \ 89 SWAP1(X[l],X[r],t1); \ 90 l++; \ 91 r--; \ 92 } \ 93 } while(0) 94 95 #define TwoWayPartition2(X,Y,pivot,t1,t2,lo,hi,l,r) \ 96 do { \ 97 l = lo; \ 98 r = hi; \ 99 while(1) { \ 100 while (X[l] < pivot) l++; \ 101 while (X[r] > pivot) r--; \ 102 if (l >= r) {r++; break;} \ 103 SWAP2(X[l],X[r],Y[l],Y[r],t1,t2); \ 104 l++; \ 105 r--; \ 106 } \ 107 } while(0) 108 109 #define TwoWayPartition3(X,Y,Z,pivot,t1,t2,t3,lo,hi,l,r) \ 110 do { \ 111 l = lo; \ 112 r = hi; \ 113 while(1) { \ 114 while (X[l] < pivot) l++; \ 115 while (X[r] > pivot) r--; \ 116 if (l >= r) {r++; break;} \ 117 SWAP3(X[l],X[r],Y[l],Y[r],Z[l],Z[r],t1,t2,t3); \ 118 l++; \ 119 r--; \ 120 } \ 121 } while(0) 122 123 /* Templates for similar functions used below */ 124 #define QuickSort1(FuncName,X,n,pivot,t1,ierr) \ 125 do { \ 126 PetscInt i,j,p,l,r,hi=n-1; \ 127 if (n < 8) { \ 128 for (i=0; i<n; i++) { \ 129 pivot = X[i]; \ 130 for (j=i+1; j<n; j++) { \ 131 if (pivot > X[j]) { \ 132 SWAP1(X[i],X[j],t1); \ 133 pivot = X[i]; \ 134 } \ 135 } \ 136 } \ 137 } else { \ 138 p = MEDIAN(X,hi); \ 139 pivot = X[p]; \ 140 TwoWayPartition1(X,pivot,t1,0,hi,l,r); \ 141 ierr = FuncName(l,X);CHKERRQ(ierr); \ 142 ierr = FuncName(hi-r+1,X+r);CHKERRQ(ierr); \ 143 } \ 144 } while(0) 145 146 /* Templates for similar functions used below */ 147 #define QuickSortReverse1(FuncName,X,n,pivot,t1,ierr) \ 148 do { \ 149 PetscInt i,j,p,l,r,hi=n-1; \ 150 if (n < 8) { \ 151 for (i=0; i<n; i++) { \ 152 pivot = X[i]; \ 153 for (j=i+1; j<n; j++) { \ 154 if (pivot < X[j]) { \ 155 SWAP1(X[i],X[j],t1); \ 156 pivot = X[i]; \ 157 } \ 158 } \ 159 } \ 160 } else { \ 161 p = MEDIAN(X,hi); \ 162 pivot = X[p]; \ 163 TwoWayPartitionReverse1(X,pivot,t1,0,hi,l,r); \ 164 ierr = FuncName(l,X);CHKERRQ(ierr); \ 165 ierr = FuncName(hi-r+1,X+r);CHKERRQ(ierr); \ 166 } \ 167 } while(0) 168 169 #define QuickSort2(FuncName,X,Y,n,pivot,t1,t2,ierr) \ 170 do { \ 171 PetscInt i,j,p,l,r,hi=n-1; \ 172 if (n < 8) { \ 173 for (i=0; i<n; i++) { \ 174 pivot = X[i]; \ 175 for (j=i+1; j<n; j++) { \ 176 if (pivot > X[j]) { \ 177 SWAP2(X[i],X[j],Y[i],Y[j],t1,t2); \ 178 pivot = X[i]; \ 179 } \ 180 } \ 181 } \ 182 } else { \ 183 p = MEDIAN(X,hi); \ 184 pivot = X[p]; \ 185 TwoWayPartition2(X,Y,pivot,t1,t2,0,hi,l,r); \ 186 ierr = FuncName(l,X,Y);CHKERRQ(ierr); \ 187 ierr = FuncName(hi-r+1,X+r,Y+r);CHKERRQ(ierr); \ 188 } \ 189 } while(0) 190 191 #define QuickSort3(FuncName,X,Y,Z,n,pivot,t1,t2,t3,ierr) \ 192 do { \ 193 PetscInt i,j,p,l,r,hi=n-1; \ 194 if (n < 8) { \ 195 for (i=0; i<n; i++) { \ 196 pivot = X[i]; \ 197 for (j=i+1; j<n; j++) { \ 198 if (pivot > X[j]) { \ 199 SWAP3(X[i],X[j],Y[i],Y[j],Z[i],Z[j],t1,t2,t3); \ 200 pivot = X[i]; \ 201 } \ 202 } \ 203 } \ 204 } else { \ 205 p = MEDIAN(X,hi); \ 206 pivot = X[p]; \ 207 TwoWayPartition3(X,Y,Z,pivot,t1,t2,t3,0,hi,l,r); \ 208 ierr = FuncName(l,X,Y,Z);CHKERRQ(ierr); \ 209 ierr = FuncName(hi-r+1,X+r,Y+r,Z+r);CHKERRQ(ierr); \ 210 } \ 211 } while(0) 212 213 /*@ 214 PetscSortInt - Sorts an array of integers in place in increasing order. 215 216 Not Collective 217 218 Input Parameters: 219 + n - number of values 220 - X - array of integers 221 222 Level: intermediate 223 224 .seealso: PetscSortReal(), PetscSortIntWithPermutation() 225 @*/ 226 PetscErrorCode PetscSortInt(PetscInt n,PetscInt X[]) 227 { 228 PetscErrorCode ierr; 229 PetscInt pivot,t1; 230 231 PetscFunctionBegin; 232 QuickSort1(PetscSortInt,X,n,pivot,t1,ierr); 233 PetscFunctionReturn(0); 234 } 235 236 /*@ 237 PetscSortReverseInt - Sorts an array of integers in place in decreasing order. 238 239 Not Collective 240 241 Input Parameters: 242 + n - number of values 243 - X - array of integers 244 245 Level: intermediate 246 247 .seealso: PetscSortInt(), PetscSortIntWithPermutation() 248 @*/ 249 PetscErrorCode PetscSortReverseInt(PetscInt n,PetscInt X[]) 250 { 251 PetscErrorCode ierr; 252 PetscInt pivot,t1; 253 254 PetscFunctionBegin; 255 QuickSortReverse1(PetscSortReverseInt,X,n,pivot,t1,ierr); 256 PetscFunctionReturn(0); 257 } 258 259 /*@ 260 PetscSortedRemoveDupsInt - Removes all duplicate entries of a sorted input array 261 262 Not Collective 263 264 Input Parameters: 265 + n - number of values 266 - X - sorted array of integers 267 268 Output Parameter: 269 . n - number of non-redundant values 270 271 Level: intermediate 272 273 .seealso: PetscSortInt() 274 @*/ 275 PetscErrorCode PetscSortedRemoveDupsInt(PetscInt *n,PetscInt X[]) 276 { 277 PetscInt i,s = 0,N = *n, b = 0; 278 279 PetscFunctionBegin; 280 for (i=0; i<N-1; i++) { 281 if (PetscUnlikely(X[b+s+1] < X[b])) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Input array is not sorted"); 282 if (X[b+s+1] != X[b]) { 283 X[b+1] = X[b+s+1]; b++; 284 } else s++; 285 } 286 *n = N - s; 287 PetscFunctionReturn(0); 288 } 289 290 /*@ 291 PetscSortRemoveDupsInt - Sorts an array of integers in place in increasing order removes all duplicate entries 292 293 Not Collective 294 295 Input Parameters: 296 + n - number of values 297 - X - array of integers 298 299 Output Parameter: 300 . n - number of non-redundant values 301 302 Level: intermediate 303 304 .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt(), PetscSortedRemoveDupsInt() 305 @*/ 306 PetscErrorCode PetscSortRemoveDupsInt(PetscInt *n,PetscInt X[]) 307 { 308 PetscErrorCode ierr; 309 310 PetscFunctionBegin; 311 ierr = PetscSortInt(*n,X);CHKERRQ(ierr); 312 ierr = PetscSortedRemoveDupsInt(n,X);CHKERRQ(ierr); 313 PetscFunctionReturn(0); 314 } 315 316 /*@ 317 PetscFindInt - Finds integer in a sorted array of integers 318 319 Not Collective 320 321 Input Parameters: 322 + key - the integer to locate 323 . n - number of values in the array 324 - X - array of integers 325 326 Output Parameter: 327 . loc - the location if found, otherwise -(slot+1) where slot is the place the value would go 328 329 Level: intermediate 330 331 .seealso: PetscSortInt(), PetscSortIntWithArray(), PetscSortRemoveDupsInt() 332 @*/ 333 PetscErrorCode PetscFindInt(PetscInt key, PetscInt n, const PetscInt X[], PetscInt *loc) 334 { 335 PetscInt lo = 0,hi = n; 336 337 PetscFunctionBegin; 338 PetscValidPointer(loc,4); 339 if (!n) {*loc = -1; PetscFunctionReturn(0);} 340 PetscValidPointer(X,3); 341 while (hi - lo > 1) { 342 PetscInt mid = lo + (hi - lo)/2; 343 if (key < X[mid]) hi = mid; 344 else lo = mid; 345 } 346 *loc = key == X[lo] ? lo : -(lo + (key > X[lo]) + 1); 347 PetscFunctionReturn(0); 348 } 349 350 /*@ 351 PetscCheckDupsInt - Checks if an integer array has duplicates 352 353 Not Collective 354 355 Input Parameters: 356 + n - number of values in the array 357 - X - array of integers 358 359 360 Output Parameter: 361 . dups - True if the array has dups, otherwise false 362 363 Level: intermediate 364 365 .seealso: PetscSortRemoveDupsInt() 366 @*/ 367 PetscErrorCode PetscCheckDupsInt(PetscInt n,const PetscInt X[],PetscBool *dups) 368 { 369 PetscErrorCode ierr; 370 PetscInt i; 371 PetscHSetI ht; 372 PetscBool missing; 373 374 PetscFunctionBegin; 375 PetscValidPointer(dups,3); 376 *dups = PETSC_FALSE; 377 if (n > 1) { 378 ierr = PetscHSetICreate(&ht);CHKERRQ(ierr); 379 ierr = PetscHSetIResize(ht,n);CHKERRQ(ierr); 380 for (i=0; i<n; i++) { 381 ierr = PetscHSetIQueryAdd(ht,X[i],&missing);CHKERRQ(ierr); 382 if(!missing) {*dups = PETSC_TRUE; break;} 383 } 384 ierr = PetscHSetIDestroy(&ht);CHKERRQ(ierr); 385 } 386 PetscFunctionReturn(0); 387 } 388 389 /*@ 390 PetscFindMPIInt - Finds MPI integer in a sorted array of integers 391 392 Not Collective 393 394 Input Parameters: 395 + key - the integer to locate 396 . n - number of values in the array 397 - X - array of integers 398 399 Output Parameter: 400 . loc - the location if found, otherwise -(slot+1) where slot is the place the value would go 401 402 Level: intermediate 403 404 .seealso: PetscSortInt(), PetscSortIntWithArray(), PetscSortRemoveDupsInt() 405 @*/ 406 PetscErrorCode PetscFindMPIInt(PetscMPIInt key, PetscInt n, const PetscMPIInt X[], PetscInt *loc) 407 { 408 PetscInt lo = 0,hi = n; 409 410 PetscFunctionBegin; 411 PetscValidPointer(loc,4); 412 if (!n) {*loc = -1; PetscFunctionReturn(0);} 413 PetscValidPointer(X,3); 414 while (hi - lo > 1) { 415 PetscInt mid = lo + (hi - lo)/2; 416 if (key < X[mid]) hi = mid; 417 else lo = mid; 418 } 419 *loc = key == X[lo] ? lo : -(lo + (key > X[lo]) + 1); 420 PetscFunctionReturn(0); 421 } 422 423 /*@ 424 PetscSortIntWithArray - Sorts an array of integers in place in increasing order; 425 changes a second array to match the sorted first array. 426 427 Not Collective 428 429 Input Parameters: 430 + n - number of values 431 . X - array of integers 432 - Y - second array of integers 433 434 Level: intermediate 435 436 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt() 437 @*/ 438 PetscErrorCode PetscSortIntWithArray(PetscInt n,PetscInt X[],PetscInt Y[]) 439 { 440 PetscErrorCode ierr; 441 PetscInt pivot,t1,t2; 442 443 PetscFunctionBegin; 444 QuickSort2(PetscSortIntWithArray,X,Y,n,pivot,t1,t2,ierr); 445 PetscFunctionReturn(0); 446 } 447 448 /*@ 449 PetscSortIntWithArrayPair - Sorts an array of integers in place in increasing order; 450 changes a pair of integer arrays to match the sorted first array. 451 452 Not Collective 453 454 Input Parameters: 455 + n - number of values 456 . X - array of integers 457 . Y - second array of integers (first array of the pair) 458 - Z - third array of integers (second array of the pair) 459 460 Level: intermediate 461 462 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortIntWithArray() 463 @*/ 464 PetscErrorCode PetscSortIntWithArrayPair(PetscInt n,PetscInt X[],PetscInt Y[],PetscInt Z[]) 465 { 466 PetscErrorCode ierr; 467 PetscInt pivot,t1,t2,t3; 468 469 PetscFunctionBegin; 470 QuickSort3(PetscSortIntWithArrayPair,X,Y,Z,n,pivot,t1,t2,t3,ierr); 471 PetscFunctionReturn(0); 472 } 473 474 /*@ 475 PetscSortMPIInt - Sorts an array of MPI integers in place in increasing order. 476 477 Not Collective 478 479 Input Parameters: 480 + n - number of values 481 - X - array of integers 482 483 Level: intermediate 484 485 .seealso: PetscSortReal(), PetscSortIntWithPermutation() 486 @*/ 487 PetscErrorCode PetscSortMPIInt(PetscInt n,PetscMPIInt X[]) 488 { 489 PetscErrorCode ierr; 490 PetscMPIInt pivot,t1; 491 492 PetscFunctionBegin; 493 QuickSort1(PetscSortMPIInt,X,n,pivot,t1,ierr); 494 PetscFunctionReturn(0); 495 } 496 497 /*@ 498 PetscSortRemoveDupsMPIInt - Sorts an array of MPI integers in place in increasing order removes all duplicate entries 499 500 Not Collective 501 502 Input Parameters: 503 + n - number of values 504 - X - array of integers 505 506 Output Parameter: 507 . n - number of non-redundant values 508 509 Level: intermediate 510 511 .seealso: PetscSortReal(), PetscSortIntWithPermutation(), PetscSortInt() 512 @*/ 513 PetscErrorCode PetscSortRemoveDupsMPIInt(PetscInt *n,PetscMPIInt X[]) 514 { 515 PetscErrorCode ierr; 516 PetscInt i,s = 0,N = *n, b = 0; 517 518 PetscFunctionBegin; 519 ierr = PetscSortMPIInt(N,X);CHKERRQ(ierr); 520 for (i=0; i<N-1; i++) { 521 if (X[b+s+1] != X[b]) { 522 X[b+1] = X[b+s+1]; b++; 523 } else s++; 524 } 525 *n = N - s; 526 PetscFunctionReturn(0); 527 } 528 529 /*@ 530 PetscSortMPIIntWithArray - Sorts an array of integers in place in increasing order; 531 changes a second array to match the sorted first array. 532 533 Not Collective 534 535 Input Parameters: 536 + n - number of values 537 . X - array of integers 538 - Y - second array of integers 539 540 Level: intermediate 541 542 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt() 543 @*/ 544 PetscErrorCode PetscSortMPIIntWithArray(PetscMPIInt n,PetscMPIInt X[],PetscMPIInt Y[]) 545 { 546 PetscErrorCode ierr; 547 PetscMPIInt pivot,t1,t2; 548 549 PetscFunctionBegin; 550 QuickSort2(PetscSortMPIIntWithArray,X,Y,n,pivot,t1,t2,ierr); 551 PetscFunctionReturn(0); 552 } 553 554 /*@ 555 PetscSortMPIIntWithIntArray - Sorts an array of MPI integers in place in increasing order; 556 changes a second array of Petsc intergers to match the sorted first array. 557 558 Not Collective 559 560 Input Parameters: 561 + n - number of values 562 . X - array of MPI integers 563 - Y - second array of Petsc integers 564 565 Level: intermediate 566 567 Notes: this routine is useful when one needs to sort MPI ranks with other integer arrays. 568 569 .seealso: PetscSortMPIIntWithArray() 570 @*/ 571 PetscErrorCode PetscSortMPIIntWithIntArray(PetscMPIInt n,PetscMPIInt X[],PetscInt Y[]) 572 { 573 PetscErrorCode ierr; 574 PetscMPIInt pivot,t1; 575 PetscInt t2; 576 577 PetscFunctionBegin; 578 QuickSort2(PetscSortMPIIntWithIntArray,X,Y,n,pivot,t1,t2,ierr); 579 PetscFunctionReturn(0); 580 } 581 582 /*@ 583 PetscSortIntWithScalarArray - Sorts an array of integers in place in increasing order; 584 changes a second SCALAR array to match the sorted first INTEGER array. 585 586 Not Collective 587 588 Input Parameters: 589 + n - number of values 590 . X - array of integers 591 - Y - second array of scalars 592 593 Level: intermediate 594 595 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray() 596 @*/ 597 PetscErrorCode PetscSortIntWithScalarArray(PetscInt n,PetscInt X[],PetscScalar Y[]) 598 { 599 PetscErrorCode ierr; 600 PetscInt pivot,t1; 601 PetscScalar t2; 602 603 PetscFunctionBegin; 604 QuickSort2(PetscSortIntWithScalarArray,X,Y,n,pivot,t1,t2,ierr); 605 PetscFunctionReturn(0); 606 } 607 608 /*@C 609 PetscSortIntWithDataArray - Sorts an array of integers in place in increasing order; 610 changes a second array to match the sorted first INTEGER array. Unlike other sort routines, the user must 611 provide workspace (the size of an element in the data array) to use when sorting. 612 613 Not Collective 614 615 Input Parameters: 616 + n - number of values 617 . X - array of integers 618 . Y - second array of data 619 . size - sizeof elements in the data array in bytes 620 - t2 - workspace of "size" bytes used when sorting 621 622 Level: intermediate 623 624 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray() 625 @*/ 626 PetscErrorCode PetscSortIntWithDataArray(PetscInt n,PetscInt X[],void *Y,size_t size,void *t2) 627 { 628 PetscErrorCode ierr; 629 char *YY = (char*)Y; 630 PetscInt i,j,p,t1,pivot,hi=n-1,l,r; 631 632 PetscFunctionBegin; 633 if (n<8) { 634 for (i=0; i<n; i++) { 635 pivot = X[i]; 636 for (j=i+1; j<n; j++) { 637 if (pivot > X[j]) { 638 SWAP2Data(X[i],X[j],YY+size*i,YY+size*j,t1,t2,size); 639 pivot = X[i]; 640 } 641 } 642 } 643 } else { 644 /* Two way partition */ 645 p = MEDIAN(X,hi); 646 pivot = X[p]; 647 l = 0; 648 r = hi; 649 while(1) { 650 while (X[l] < pivot) l++; 651 while (X[r] > pivot) r--; 652 if (l >= r) {r++; break;} 653 SWAP2Data(X[l],X[r],YY+size*l,YY+size*r,t1,t2,size); 654 l++; 655 r--; 656 } 657 ierr = PetscSortIntWithDataArray(l,X,Y,size,t2);CHKERRQ(ierr); 658 ierr = PetscSortIntWithDataArray(hi-r+1,X+r,YY+size*r,size,t2);CHKERRQ(ierr); 659 } 660 PetscFunctionReturn(0); 661 } 662 663 /*@ 664 PetscMergeIntArray - Merges two SORTED integer arrays, removes duplicate elements. 665 666 Not Collective 667 668 Input Parameters: 669 + an - number of values in the first array 670 . aI - first sorted array of integers 671 . bn - number of values in the second array 672 - bI - second array of integers 673 674 Output Parameters: 675 + n - number of values in the merged array 676 - L - merged sorted array, this is allocated if an array is not provided 677 678 Level: intermediate 679 680 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray() 681 @*/ 682 PetscErrorCode PetscMergeIntArray(PetscInt an,const PetscInt aI[], PetscInt bn, const PetscInt bI[], PetscInt *n, PetscInt **L) 683 { 684 PetscErrorCode ierr; 685 PetscInt *L_ = *L, ak, bk, k; 686 687 if (!L_) { 688 ierr = PetscMalloc1(an+bn, L);CHKERRQ(ierr); 689 L_ = *L; 690 } 691 k = ak = bk = 0; 692 while (ak < an && bk < bn) { 693 if (aI[ak] == bI[bk]) { 694 L_[k] = aI[ak]; 695 ++ak; 696 ++bk; 697 ++k; 698 } else if (aI[ak] < bI[bk]) { 699 L_[k] = aI[ak]; 700 ++ak; 701 ++k; 702 } else { 703 L_[k] = bI[bk]; 704 ++bk; 705 ++k; 706 } 707 } 708 if (ak < an) { 709 ierr = PetscArraycpy(L_+k,aI+ak,an-ak);CHKERRQ(ierr); 710 k += (an-ak); 711 } 712 if (bk < bn) { 713 ierr = PetscArraycpy(L_+k,bI+bk,bn-bk);CHKERRQ(ierr); 714 k += (bn-bk); 715 } 716 *n = k; 717 PetscFunctionReturn(0); 718 } 719 720 /*@ 721 PetscMergeIntArrayPair - Merges two SORTED integer arrays that share NO common values along with an additional array of integers. 722 The additional arrays are the same length as sorted arrays and are merged 723 in the order determined by the merging of the sorted pair. 724 725 Not Collective 726 727 Input Parameters: 728 + an - number of values in the first array 729 . aI - first sorted array of integers 730 . aJ - first additional array of integers 731 . bn - number of values in the second array 732 . bI - second array of integers 733 - bJ - second additional of integers 734 735 Output Parameters: 736 + n - number of values in the merged array (== an + bn) 737 . L - merged sorted array 738 - J - merged additional array 739 740 Notes: 741 if L or J point to non-null arrays then this routine will assume they are of the approproate size and use them, otherwise this routine will allocate space for them 742 Level: intermediate 743 744 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray() 745 @*/ 746 PetscErrorCode PetscMergeIntArrayPair(PetscInt an,const PetscInt aI[], const PetscInt aJ[], PetscInt bn, const PetscInt bI[], const PetscInt bJ[], PetscInt *n, PetscInt **L, PetscInt **J) 747 { 748 PetscErrorCode ierr; 749 PetscInt n_, *L_, *J_, ak, bk, k; 750 751 PetscFunctionBegin; 752 PetscValidIntPointer(L,8); 753 PetscValidIntPointer(J,9); 754 n_ = an + bn; 755 *n = n_; 756 if (!*L) { 757 ierr = PetscMalloc1(n_, L);CHKERRQ(ierr); 758 } 759 L_ = *L; 760 if (!*J) { 761 ierr = PetscMalloc1(n_, J);CHKERRQ(ierr); 762 } 763 J_ = *J; 764 k = ak = bk = 0; 765 while (ak < an && bk < bn) { 766 if (aI[ak] <= bI[bk]) { 767 L_[k] = aI[ak]; 768 J_[k] = aJ[ak]; 769 ++ak; 770 ++k; 771 } else { 772 L_[k] = bI[bk]; 773 J_[k] = bJ[bk]; 774 ++bk; 775 ++k; 776 } 777 } 778 if (ak < an) { 779 ierr = PetscArraycpy(L_+k,aI+ak,an-ak);CHKERRQ(ierr); 780 ierr = PetscArraycpy(J_+k,aJ+ak,an-ak);CHKERRQ(ierr); 781 k += (an-ak); 782 } 783 if (bk < bn) { 784 ierr = PetscArraycpy(L_+k,bI+bk,bn-bk);CHKERRQ(ierr); 785 ierr = PetscArraycpy(J_+k,bJ+bk,bn-bk);CHKERRQ(ierr); 786 } 787 PetscFunctionReturn(0); 788 } 789 790 /*@ 791 PetscMergeMPIIntArray - Merges two SORTED integer arrays. 792 793 Not Collective 794 795 Input Parameters: 796 + an - number of values in the first array 797 . aI - first sorted array of integers 798 . bn - number of values in the second array 799 - bI - second array of integers 800 801 Output Parameters: 802 + n - number of values in the merged array (<= an + bn) 803 - L - merged sorted array, allocated if address of NULL pointer is passed 804 805 Level: intermediate 806 807 .seealso: PetscSortReal(), PetscSortIntPermutation(), PetscSortInt(), PetscSortIntWithArray() 808 @*/ 809 PetscErrorCode PetscMergeMPIIntArray(PetscInt an,const PetscMPIInt aI[],PetscInt bn,const PetscMPIInt bI[],PetscInt *n,PetscMPIInt **L) 810 { 811 PetscErrorCode ierr; 812 PetscInt ai,bi,k; 813 814 PetscFunctionBegin; 815 if (!*L) {ierr = PetscMalloc1((an+bn),L);CHKERRQ(ierr);} 816 for (ai=0,bi=0,k=0; ai<an || bi<bn; ) { 817 PetscInt t = -1; 818 for ( ; ai<an && (!bn || aI[ai] <= bI[bi]); ai++) (*L)[k++] = t = aI[ai]; 819 for ( ; bi<bn && bI[bi] == t; bi++); 820 for ( ; bi<bn && (!an || bI[bi] <= aI[ai]); bi++) (*L)[k++] = t = bI[bi]; 821 for ( ; ai<an && aI[ai] == t; ai++); 822 } 823 *n = k; 824 PetscFunctionReturn(0); 825 } 826 827 /*@C 828 PetscProcessTree - Prepares tree data to be displayed graphically 829 830 Not Collective 831 832 Input Parameters: 833 + n - number of values 834 . mask - indicates those entries in the tree, location 0 is always masked 835 - parentid - indicates the parent of each entry 836 837 Output Parameters: 838 + Nlevels - the number of levels 839 . Level - for each node tells its level 840 . Levelcnts - the number of nodes on each level 841 . Idbylevel - a list of ids on each of the levels, first level followed by second etc 842 - Column - for each id tells its column index 843 844 Level: developer 845 846 Notes: 847 This code is not currently used 848 849 .seealso: PetscSortReal(), PetscSortIntWithPermutation() 850 @*/ 851 PetscErrorCode PetscProcessTree(PetscInt n,const PetscBool mask[],const PetscInt parentid[],PetscInt *Nlevels,PetscInt **Level,PetscInt **Levelcnt,PetscInt **Idbylevel,PetscInt **Column) 852 { 853 PetscInt i,j,cnt,nmask = 0,nlevels = 0,*level,*levelcnt,levelmax = 0,*workid,*workparentid,tcnt = 0,*idbylevel,*column; 854 PetscErrorCode ierr; 855 PetscBool done = PETSC_FALSE; 856 857 PetscFunctionBegin; 858 if (!mask[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Mask of 0th location must be set"); 859 for (i=0; i<n; i++) { 860 if (mask[i]) continue; 861 if (parentid[i] == i) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Node labeled as own parent"); 862 if (parentid[i] && mask[parentid[i]]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Parent is masked"); 863 } 864 865 for (i=0; i<n; i++) { 866 if (!mask[i]) nmask++; 867 } 868 869 /* determine the level in the tree of each node */ 870 ierr = PetscCalloc1(n,&level);CHKERRQ(ierr); 871 872 level[0] = 1; 873 while (!done) { 874 done = PETSC_TRUE; 875 for (i=0; i<n; i++) { 876 if (mask[i]) continue; 877 if (!level[i] && level[parentid[i]]) level[i] = level[parentid[i]] + 1; 878 else if (!level[i]) done = PETSC_FALSE; 879 } 880 } 881 for (i=0; i<n; i++) { 882 level[i]--; 883 nlevels = PetscMax(nlevels,level[i]); 884 } 885 886 /* count the number of nodes on each level and its max */ 887 ierr = PetscCalloc1(nlevels,&levelcnt);CHKERRQ(ierr); 888 for (i=0; i<n; i++) { 889 if (mask[i]) continue; 890 levelcnt[level[i]-1]++; 891 } 892 for (i=0; i<nlevels;i++) levelmax = PetscMax(levelmax,levelcnt[i]); 893 894 /* for each level sort the ids by the parent id */ 895 ierr = PetscMalloc2(levelmax,&workid,levelmax,&workparentid);CHKERRQ(ierr); 896 ierr = PetscMalloc1(nmask,&idbylevel);CHKERRQ(ierr); 897 for (j=1; j<=nlevels;j++) { 898 cnt = 0; 899 for (i=0; i<n; i++) { 900 if (mask[i]) continue; 901 if (level[i] != j) continue; 902 workid[cnt] = i; 903 workparentid[cnt++] = parentid[i]; 904 } 905 /* PetscIntView(cnt,workparentid,0); 906 PetscIntView(cnt,workid,0); 907 ierr = PetscSortIntWithArray(cnt,workparentid,workid);CHKERRQ(ierr); 908 PetscIntView(cnt,workparentid,0); 909 PetscIntView(cnt,workid,0);*/ 910 ierr = PetscArraycpy(idbylevel+tcnt,workid,cnt);CHKERRQ(ierr); 911 tcnt += cnt; 912 } 913 if (tcnt != nmask) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent count of unmasked nodes"); 914 ierr = PetscFree2(workid,workparentid);CHKERRQ(ierr); 915 916 /* for each node list its column */ 917 ierr = PetscMalloc1(n,&column);CHKERRQ(ierr); 918 cnt = 0; 919 for (j=0; j<nlevels; j++) { 920 for (i=0; i<levelcnt[j]; i++) { 921 column[idbylevel[cnt++]] = i; 922 } 923 } 924 925 *Nlevels = nlevels; 926 *Level = level; 927 *Levelcnt = levelcnt; 928 *Idbylevel = idbylevel; 929 *Column = column; 930 PetscFunctionReturn(0); 931 } 932 933 /*@ 934 PetscParallelSortedInt - Check whether an integer array, distributed over a communicator, is globally sorted. 935 936 Collective 937 938 Input Parameters: 939 + comm - the MPI communicator 940 . n - the local number of integers 941 - keys - the local array of integers 942 943 Output Parameters: 944 . is_sorted - whether the array is globally sorted 945 946 Level: developer 947 948 .seealso: PetscParallelSortInt() 949 @*/ 950 PetscErrorCode PetscParallelSortedInt(MPI_Comm comm, PetscInt n, const PetscInt keys[], PetscBool *is_sorted) 951 { 952 PetscBool sorted; 953 PetscInt i, min, max, prevmax; 954 PetscMPIInt rank; 955 PetscErrorCode ierr; 956 957 PetscFunctionBegin; 958 sorted = PETSC_TRUE; 959 min = PETSC_MAX_INT; 960 max = PETSC_MIN_INT; 961 if (n) { 962 min = keys[0]; 963 max = keys[0]; 964 } 965 for (i = 1; i < n; i++) { 966 if (keys[i] < keys[i - 1]) break; 967 min = PetscMin(min,keys[i]); 968 max = PetscMax(max,keys[i]); 969 } 970 if (i < n) sorted = PETSC_FALSE; 971 prevmax = PETSC_MIN_INT; 972 ierr = MPI_Exscan(&max, &prevmax, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr); 973 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 974 if (!rank) prevmax = PETSC_MIN_INT; 975 if (prevmax > min) sorted = PETSC_FALSE; 976 ierr = MPI_Allreduce(&sorted, is_sorted, 1, MPIU_BOOL, MPI_LAND, comm);CHKERRQ(ierr); 977 PetscFunctionReturn(0); 978 } 979 980