1 //------------------------------------------------------------------------------ 2 // GB_sparse_masker_template: R = masker (C, M, Z) where R is sparse/hyper 3 //------------------------------------------------------------------------------ 4 5 // SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2021, All Rights Reserved. 6 // SPDX-License-Identifier: Apache-2.0 7 8 //------------------------------------------------------------------------------ 9 10 // Computes C<M>=Z or C<!M>=Z, returning the result in R, which is sparse or 11 // hypersparse. The input matrix C is not modified. Effectively, this 12 // computes R=C and then R<M>=Z or R<!M>=Z. If the C_replace descriptor is 13 // enabled, then C has already been cleared, and is an empty (but non-NULL) 14 // matrix. 15 16 // phase1: does not compute R itself, but just counts the # of entries in each 17 // vector of R. Fine tasks compute the # of entries in their slice of a 18 // single vector of R, and the results are cumsum'd. 19 20 // phase2: computes R, using the counts computed by phase1. 21 22 // C is sparse or hypersparse. M and Z can have any sparsity structure. 23 24 // ------------------------------------------ 25 // C <!M> = Z R 26 // ------------------------------------------ 27 28 // sparse sparse sparse sparse 29 // sparse bitmap sparse sparse 30 // sparse full sparse sparse 31 32 // ------------------------------------------ 33 // C <M> = Z R 34 // ------------------------------------------ 35 36 // sparse sparse sparse sparse 37 // sparse sparse bitmap sparse 38 // sparse sparse full sparse 39 // sparse bitmap sparse sparse 40 // sparse full sparse sparse 41 42 // FUTURE:: add special cases for C==Z, C==M, and Z==M aliases 43 44 //------------------------------------------------------------------------------ 45 // R(i,j) = Z(i,j) when Z is sparse or hypersparse 46 //------------------------------------------------------------------------------ 47 48 #if defined ( GB_PHASE_1_OF_2 ) 49 #define GB_COPY_Z \ 50 { \ 51 rjnz++ ; \ 52 } 53 #else 54 #define GB_COPY_Z \ 55 { \ 56 Ri [pR] = i ; \ 57 memcpy (Rx +(pR)*rsize, Zx +(pZ)*rsize, rsize) ; \ 58 pR++ ; \ 59 } 60 #endif 61 62 //------------------------------------------------------------------------------ 63 // R(i,j) = Z(i,j) when Z is bitmap or full 64 //------------------------------------------------------------------------------ 65 66 #if defined ( GB_PHASE_1_OF_2 ) 67 #define GB_COPY_Z_BITMAP_OR_FULL \ 68 { \ 69 rjnz += GBB (Zb, pZ_start + i - iZ_first) ; \ 70 } 71 #else 72 #define GB_COPY_Z_BITMAP_OR_FULL \ 73 { \ 74 int64_t pZ = pZ_start + i - iZ_first ; \ 75 if (GBB (Zb, pZ)) \ 76 { \ 77 Ri [pR] = i ; \ 78 memcpy (Rx +(pR)*rsize, Zx +(pZ)*rsize, rsize) ; \ 79 pR++ ; \ 80 } \ 81 } 82 #endif 83 84 //------------------------------------------------------------------------------ 85 // R(i,j) = C(i,j) 86 //------------------------------------------------------------------------------ 87 88 #if defined ( GB_PHASE_1_OF_2 ) 89 #define GB_COPY_C \ 90 { \ 91 rjnz++ ; \ 92 } 93 #else 94 #define GB_COPY_C \ 95 { \ 96 Ri [pR] = i ; \ 97 memcpy (Rx +(pR)*rsize, Cx +(pC)*rsize, rsize) ; \ 98 pR++ ; \ 99 } 100 #endif 101 102 //------------------------------------------------------------------------------ 103 // template for R = masker (C, M, Z) when R is sparse or hypersparse 104 //------------------------------------------------------------------------------ 105 106 { 107 108 //-------------------------------------------------------------------------- 109 // phase1: count entries in each C(:,j) 110 // phase2: compute C 111 //-------------------------------------------------------------------------- 112 113 ASSERT (C_is_sparse || C_is_hyper) ; 114 115 #pragma omp parallel for num_threads(R_nthreads) schedule(dynamic,1) 116 for (taskid = 0 ; taskid < R_ntasks ; taskid++) 117 { 118 119 //---------------------------------------------------------------------- 120 // get the task descriptor 121 //---------------------------------------------------------------------- 122 123 int64_t kfirst = TaskList [taskid].kfirst ; 124 int64_t klast = TaskList [taskid].klast ; 125 bool fine_task = (klast == -1) ; 126 int64_t len ; 127 if (fine_task) 128 { 129 // a fine task operates on a slice of a single vector 130 klast = kfirst ; 131 len = TaskList [taskid].len ; 132 } 133 else 134 { 135 // a coarse task operates on one or more whole vectors 136 len = vlen ; 137 } 138 139 //---------------------------------------------------------------------- 140 // compute all vectors in this task 141 //---------------------------------------------------------------------- 142 143 for (int64_t k = kfirst ; k <= klast ; k++) 144 { 145 146 //------------------------------------------------------------------ 147 // get j, the kth vector of R 148 //------------------------------------------------------------------ 149 150 int64_t j = GBH (Rh, k) ; 151 152 #if defined ( GB_PHASE_1_OF_2 ) 153 int64_t rjnz = 0 ; 154 #else 155 int64_t pR, pR_end ; 156 if (fine_task) 157 { 158 // A fine task computes a slice of R(:,j) 159 pR = TaskList [taskid ].pC ; 160 pR_end = TaskList [taskid+1].pC ; 161 ASSERT (Rp [k] <= pR && pR <= pR_end && pR_end <= Rp [k+1]) ; 162 } 163 else 164 { 165 // The vectors of R are never sliced for a coarse task. 166 pR = Rp [k] ; 167 pR_end = Rp [k+1] ; 168 } 169 int64_t rjnz = pR_end - pR ; 170 if (rjnz == 0) 171 { 172 continue ; 173 } 174 #endif 175 176 //------------------------------------------------------------------ 177 // get C(:,j) 178 //------------------------------------------------------------------ 179 180 int64_t pC = -1, pC_end = -1 ; 181 if (fine_task) 182 { 183 // A fine task operates on Ci,Cx [pC...pC_end-1], which is 184 // a subset of the vector C(:,j) 185 pC = TaskList [taskid].pA ; 186 pC_end = TaskList [taskid].pA_end ; 187 } 188 else 189 { 190 // A coarse task operates on the entire vector C(:,j) 191 int64_t kC = (R_to_C == NULL) ? j : R_to_C [k] ; 192 if (kC >= 0) 193 { 194 pC = Cp [kC] ; 195 pC_end = Cp [kC+1] ; 196 } 197 } 198 199 int64_t cjnz = pC_end - pC ; // nnz in C(:,j) for this slice 200 bool cdense = (cjnz == len) && (cjnz > 0) ; 201 202 #if defined ( GB_PHASE_2_OF_2 ) || defined ( GB_DEBUG ) 203 // get the first index in C(:,j) for this vector 204 int64_t iC_first = -1 ; 205 if (cjnz > 0) iC_first = Ci [pC] ; 206 #endif 207 208 #ifdef GB_DEBUG 209 int64_t iC_last = -1 ; 210 if (cjnz > 0) iC_last = Ci [pC_end-1] ; 211 #endif 212 213 //------------------------------------------------------------------ 214 // get Z(:,j) 215 //------------------------------------------------------------------ 216 217 int64_t pZ = -1, pZ_end = -1 ; 218 if (fine_task) 219 { 220 // A fine task operates on Zi,Zx [pZ...pZ_end-1], which is 221 // a subset of the vector Z(:,j) 222 pZ = TaskList [taskid].pB ; 223 pZ_end = TaskList [taskid].pB_end ; 224 } 225 else 226 { 227 // A coarse task operates on the entire vector Z(:,j) 228 int64_t kZ = (R_to_Z == NULL) ? j : R_to_Z [k] ; 229 if (kZ >= 0) 230 { 231 pZ = GBP (Zp, kZ, vlen) ; 232 pZ_end = GBP (Zp, kZ+1, vlen) ; 233 } 234 } 235 236 int64_t zjnz = pZ_end - pZ ; // nnz in Z(:,j) for this slice 237 int64_t pZ_start = pZ ; 238 bool zdense = (zjnz == len) && (zjnz > 0) ; 239 240 int64_t iZ_first = -1, iZ_last = -1 ; 241 if (zjnz > 0) 242 { 243 iZ_first = GBI (Zi, pZ, vlen) ; 244 iZ_last = GBI (Zi, pZ_end-1, vlen) ; 245 } 246 247 //------------------------------------------------------------------ 248 // get M(:,j) 249 //------------------------------------------------------------------ 250 251 int64_t pM = -1, pM_end = -1 ; 252 if (fine_task) 253 { 254 // A fine task operates on Mi,Mx [pM...pM_end-1], which is 255 // a subset of the vector M(:,j) 256 pM = TaskList [taskid].pM ; 257 pM_end = TaskList [taskid].pM_end ; 258 } 259 else 260 { 261 // A coarse task operates on the entire vector M (:,j) 262 int64_t kM = (R_to_M == NULL) ? j : R_to_M [k] ; 263 if (kM >= 0) 264 { 265 pM = GBP (Mp, kM, vlen) ; 266 pM_end = GBP (Mp, kM+1, vlen) ; 267 } 268 } 269 270 int64_t mjnz = pM_end - pM ; // nnz (M (:,j)) 271 bool mdense = (mjnz == len) && (mjnz > 0) ; 272 273 // get the first index in M(:,j) for this vector 274 int64_t iM_first = -1 ; 275 int64_t pM_first = pM ; 276 if (mjnz > 0) iM_first = GBI (Mi, pM_first, vlen) ; 277 278 //------------------------------------------------------------------ 279 // R(:,j) = masker (C (:,j), M (:,j), Z (:,j)) 280 //------------------------------------------------------------------ 281 282 if (Z_is_bitmap || Z_is_full) 283 { 284 285 //-------------------------------------------------------------- 286 // Method01: Z is bitmap or full; M is sparse or hypersparse 287 //-------------------------------------------------------------- 288 289 // ------------------------------------------ 290 // C <M> = Z R 291 // ------------------------------------------ 292 293 // sparse sparse bitmap sparse 294 // sparse sparse full sparse 295 296 // M is sparse or hypersparse, and not complemented. 297 // Otherwise, R is bitmap and not computed here, but in 298 // GB_bitmap_masker_template instead. 299 300 ASSERT (M_is_sparse || M_is_hyper) ; 301 ASSERT (!Mask_comp) ; 302 303 // 2-way merge of C(:,j) and M(:,j) and direct lookup of Z 304 305 while (pC < pC_end && pM < pM_end) 306 { 307 308 int64_t iC = Ci [pC] ; 309 int64_t iM = Mi [pM] ; 310 311 if (iC < iM) 312 { 313 // C(i,j) is present but M(i,j) is not 314 // R(i,j) = C(i,j) 315 int64_t i = iC ; 316 GB_COPY_C ; 317 pC++ ; 318 } 319 else if (iC > iM) 320 { 321 // M(i,j) is present but C(i,j) is not 322 int64_t i = iM ; 323 bool mij = GB_mcast (Mx, pM, msize) ; 324 if (mij) 325 { 326 // R(i,j) = Z(i,j) 327 GB_COPY_Z_BITMAP_OR_FULL ; 328 } 329 pM++ ; 330 } 331 else 332 { 333 // both C(i,j) and M(i,j) are present 334 int64_t i = iM ; 335 bool mij = GB_mcast (Mx, pM, msize) ; 336 if (mij) 337 { 338 // R(i,j) = Z(i,j) 339 GB_COPY_Z_BITMAP_OR_FULL ; 340 } 341 else 342 { 343 // R(i,j) = C(i,j) 344 GB_COPY_C ; 345 } 346 pC++ ; 347 pM++ ; 348 } 349 } 350 351 // if M(:,j) is exhausted ; continue scanning all of C(:,j) 352 #if defined ( GB_PHASE_1_OF_2 ) 353 rjnz += (pC_end - pC) ; 354 #else 355 for ( ; pC < pC_end ; pC++) 356 { 357 // C(i,j) is present but M(i,j) is not 358 int64_t i = Ci [pC] ; 359 GB_COPY_C ; 360 } 361 #endif 362 363 // if C(:,j) is exhausted ; continue scanning all of M(:,j) 364 for ( ; pM < pM_end ; pM++) 365 { 366 // M(i,j) is present but C(i,j) is not 367 int64_t i = Mi [pM] ; 368 bool mij = GB_mcast (Mx, pM, msize) ; 369 if (mij) 370 { 371 // R(i,j) = Z(i,j) 372 GB_COPY_Z_BITMAP_OR_FULL ; 373 } 374 } 375 376 } 377 else if (mjnz == 0) 378 { 379 380 //-------------------------------------------------------------- 381 // Z is sparse or hypersparse, M(:,j) is empty 382 //-------------------------------------------------------------- 383 384 // ------------------------------------------ 385 // C <!M> = Z R 386 // ------------------------------------------ 387 388 // sparse sparse sparse sparse 389 390 // ------------------------------------------ 391 // C <M> = Z R 392 // ------------------------------------------ 393 394 // sparse sparse sparse sparse 395 396 // Z must be sparse or hypersparse 397 ASSERT (Z_is_sparse || Z_is_hyper) ; 398 399 if (!Mask_comp) 400 { 401 402 //---------------------------------------------------------- 403 // Method02: M(:,j) is empty and not complemented 404 //---------------------------------------------------------- 405 406 // R(:,j) = C(:,j), regardless of Z(:,j) 407 #if defined ( GB_PHASE_1_OF_2 ) 408 rjnz = cjnz ; 409 #else 410 ASSERT (rjnz == cjnz) ; 411 memcpy (Ri +(pR), Ci +(pC), cjnz * sizeof (int64_t)) ; 412 memcpy (Rx +(pR)*rsize, Cx +(pC)*rsize, cjnz*rsize) ; 413 #endif 414 415 } 416 else 417 { 418 419 //---------------------------------------------------------- 420 // Method03: M(:,j) is empty and complemented 421 //---------------------------------------------------------- 422 423 // R(:,j) = Z(:,j), regardless of C(:,j) 424 #if defined ( GB_PHASE_1_OF_2 ) 425 rjnz = zjnz ; 426 #else 427 ASSERT (rjnz == zjnz) ; 428 memcpy (Ri +(pR), Zi +(pZ), zjnz * sizeof (int64_t)) ; 429 memcpy (Rx +(pR)*rsize, Zx +(pZ)*rsize, zjnz*rsize) ; 430 #endif 431 } 432 433 } 434 else if (cdense && zdense) 435 { 436 437 //-------------------------------------------------------------- 438 // Method03: C(:,j) and Z(:,j) dense: thus R(:,j) dense 439 //-------------------------------------------------------------- 440 441 // ------------------------------------------ 442 // C <!M> = Z R 443 // ------------------------------------------ 444 445 // sparse sparse sparse sparse 446 // sparse bitmap sparse sparse 447 // sparse full sparse sparse 448 449 // ------------------------------------------ 450 // C <M> = Z R 451 // ------------------------------------------ 452 453 // sparse sparse sparse sparse 454 // sparse bitmap sparse sparse 455 // sparse full sparse sparse 456 457 // Both C(:,j) and Z(:,j) are dense (that is, all entries 458 // present), but both C and Z are stored in a sparse or 459 // hypersparse sparsity structure. M has any sparsity. 460 461 ASSERT (Z_is_sparse || Z_is_hyper) ; 462 463 ASSERT (cjnz == zjnz) ; 464 ASSERT (iC_first == iZ_first) ; 465 ASSERT (iC_last == iZ_last ) ; 466 #if defined ( GB_PHASE_1_OF_2 ) 467 rjnz = cjnz ; 468 #else 469 ASSERT (rjnz == cjnz) ; 470 for (int64_t p = 0 ; p < cjnz ; p++) 471 { 472 int64_t i = p + iC_first ; 473 Ri [pR + p] = i ; 474 int64_t iM = (pM < pM_end) ? GBI (Mi, pM, vlen) : INT64_MAX; 475 bool mij = false ; 476 if (i == iM) 477 { 478 mij = GBB (Mb, pM) && GB_mcast (Mx, pM, msize) ; 479 pM++ ; 480 } 481 if (Mask_comp) mij = !mij ; 482 if (mij) 483 { 484 // R(i,j) = Z (i,j) 485 memcpy (Rx +(pR+p)*rsize, Zx +(pZ+p)*rsize, rsize) ; 486 } 487 else 488 { 489 // R(i,j) = C (i,j) 490 memcpy (Rx +(pR+p)*rsize, Cx +(pC+p)*rsize, rsize) ; 491 } 492 } 493 #endif 494 495 } 496 else 497 { 498 499 //-------------------------------------------------------------- 500 // Method04: 2-way merge of C(:,j) and Z(:,j) 501 //-------------------------------------------------------------- 502 503 // Z is sparse or hypersparse; M has any sparsity structure 504 ASSERT (Z_is_sparse || Z_is_hyper) ; 505 506 //-------------------------------------------------------------- 507 // Z is sparse or hypersparse, M has any sparsity 508 //-------------------------------------------------------------- 509 510 // ------------------------------------------ 511 // C <!M> = Z R 512 // ------------------------------------------ 513 514 // sparse sparse sparse sparse 515 // sparse bitmap sparse sparse 516 // sparse full sparse sparse 517 518 // ------------------------------------------ 519 // C <M> = Z R 520 // ------------------------------------------ 521 522 // sparse sparse sparse sparse 523 // sparse bitmap sparse sparse 524 // sparse full sparse sparse 525 526 while (pC < pC_end && pZ < pZ_end) 527 { 528 529 //---------------------------------------------------------- 530 // get the next i for R(:,j) 531 //---------------------------------------------------------- 532 533 int64_t iC = Ci [pC] ; 534 int64_t iZ = Zi [pZ] ; 535 int64_t i = GB_IMIN (iC, iZ) ; 536 537 //---------------------------------------------------------- 538 // get M(i,j) 539 //---------------------------------------------------------- 540 541 bool mij = false ; 542 543 if (mdense) 544 { 545 546 //------------------------------------------------------ 547 // Method04a: M(:,j) is dense 548 //------------------------------------------------------ 549 550 // mask is dense, lookup M(i,j) 551 // iM_first == Mi [pM_first] 552 // iM_first + delta == Mi [pM_first + delta] 553 // let i = iM_first + delta 554 // let pM = pM_first + delta 555 // then delta = i - iM_first 556 pM = pM_first + (i - iM_first) ; 557 ASSERT (i == GBI (Mi, pM, vlen)) ; 558 mij = GBB (Mb, pM) && GB_mcast (Mx, pM, msize) ; 559 // increment pM for the wrapup phase below 560 pM++ ; 561 562 } 563 else 564 { 565 566 //------------------------------------------------------ 567 // Method04b: M(:,j) is sparse 568 //------------------------------------------------------ 569 570 // Use GB_SPLIT_BINARY_SEARCH so that pM can be used in 571 // the for loop with index pM in the wrapup phase. 572 ASSERT (M_is_sparse || M_is_hyper) ; 573 int64_t pright = pM_end - 1 ; 574 bool found ; 575 GB_SPLIT_BINARY_SEARCH (i, Mi, pM, pright, found) ; 576 if (found) 577 { 578 ASSERT (i == Mi [pM]) ; 579 mij = GB_mcast (Mx, pM, msize) ; 580 // increment pM for the wrapup phase below 581 pM++ ; 582 } 583 } 584 585 if (Mask_comp) mij = !mij ; 586 587 //---------------------------------------------------------- 588 // R(i,j) = C(i,j) or Z(i,j) 589 //---------------------------------------------------------- 590 591 if (iC < iZ) 592 { 593 // C(i,j) is present but Z(i,j) is not 594 if (!mij) GB_COPY_C ; 595 pC++ ; 596 } 597 else if (iC > iZ) 598 { 599 // Z(i,j) is present but C(i,j) is not 600 if (mij) GB_COPY_Z ; 601 pZ++ ; 602 } 603 else 604 { 605 // both C(i,j) and Z(i,j) are present 606 int64_t i = iC ; 607 if (mij) 608 { 609 GB_COPY_Z ; 610 } 611 else 612 { 613 GB_COPY_C ; 614 } 615 pC++ ; 616 pZ++ ; 617 } 618 } 619 620 //-------------------------------------------------------------- 621 // Method04: wrapup: C or Z are exhausted, or initially empty 622 //-------------------------------------------------------------- 623 624 cjnz = pC_end - pC ; // nnz (C(:,j)) remaining 625 zjnz = pZ_end - pZ ; // nnz (Z(:,j)) remaining 626 mjnz = pM_end - pM ; // nnz (M(:,j)) remaining 627 628 if (cjnz == 0) 629 { 630 631 //---------------------------------------------------------- 632 // C(:,j) is empty 633 //---------------------------------------------------------- 634 635 if (!Mask_comp) 636 { 637 638 //------------------------------------------------------ 639 // mask is not complemented 640 //------------------------------------------------------ 641 642 if (mdense) 643 { 644 645 //-------------------------------------------------- 646 // Method04c: M(:,j) is dense 647 //-------------------------------------------------- 648 649 for ( ; pZ < pZ_end ; pZ++) 650 { 651 int64_t i = Zi [pZ] ; 652 // mask is dense, lookup M(i,j) 653 pM = pM_first + (i - iM_first) ; 654 ASSERT (i == GBI (Mi, pM, vlen)) ; 655 bool mij = GBB (Mb, pM) && 656 GB_mcast (Mx, pM, msize) ; 657 if (mij) GB_COPY_Z ; 658 } 659 660 } 661 else if (zjnz > 32 * mjnz) 662 { 663 664 //-------------------------------------------------- 665 // Method04d: Z(:,j) is much denser than M(:,j) 666 //-------------------------------------------------- 667 668 // This loop requires pM to start at the first 669 // entry in M(:,j) that has not yet been handled. 670 671 ASSERT (M_is_sparse || M_is_hyper) ; 672 for ( ; pM < pM_end ; pM++) 673 { 674 if (GB_mcast (Mx, pM, msize)) 675 { 676 int64_t i = Mi [pM] ; 677 int64_t pright = pZ_end - 1 ; 678 bool found ; 679 GB_BINARY_SEARCH (i, Zi, pZ, pright, found); 680 if (found) GB_COPY_Z ; 681 } 682 } 683 684 } 685 else if (mjnz > 32 * zjnz) 686 { 687 688 //-------------------------------------------------- 689 // Method04e: M(:,j) is much denser than Z(:,j) 690 //-------------------------------------------------- 691 692 ASSERT (M_is_sparse || M_is_hyper) ; 693 for ( ; pZ < pZ_end ; pZ++) 694 { 695 int64_t i = Zi [pZ] ; 696 bool mij = false ; 697 int64_t pright = pM_end - 1 ; 698 bool found ; 699 GB_BINARY_SEARCH (i, Mi, pM, pright,found) ; 700 if (found) mij = GB_mcast (Mx, pM, msize) ; 701 if (mij) GB_COPY_Z ; 702 } 703 704 } 705 else 706 { 707 708 //-------------------------------------------------- 709 // Method04f: M(:,j) and Z(:,j) about same # entries 710 //-------------------------------------------------- 711 712 ASSERT (M_is_sparse || M_is_hyper) ; 713 while (pM < pM_end && pZ < pZ_end) 714 { 715 int64_t iM = Mi [pM] ; 716 int64_t i = Zi [pZ] ; 717 if (iM < i) 718 { 719 // M(i,j) exists but not Z(i,j) 720 pM++ ; 721 } 722 else if (i < iM) 723 { 724 // Z(i,j) exists but not M(i,j) 725 pZ++ ; 726 } 727 else 728 { 729 // both M(i,j) and Z(i,j) exist 730 if (GB_mcast (Mx, pM, msize)) GB_COPY_Z ; 731 pM++ ; 732 pZ++ ; 733 } 734 } 735 } 736 737 } 738 else 739 { 740 741 //------------------------------------------------------ 742 // complemented mask, and C(:,j) empty 743 //------------------------------------------------------ 744 745 if (mdense) 746 { 747 748 //-------------------------------------------------- 749 // Method04g: M(:,j) is dense 750 //-------------------------------------------------- 751 752 for ( ; pZ < pZ_end ; pZ++) 753 { 754 int64_t i = Zi [pZ] ; 755 // mask is dense, lookup M(i,j) 756 pM = pM_first + (i - iM_first) ; 757 ASSERT (i == GBI (Mi, pM, vlen)) ; 758 bool mij = GBB (Mb, pM) && 759 GB_mcast (Mx, pM, msize) ; 760 if (!mij) GB_COPY_Z ; // mask is complemented 761 } 762 } 763 else 764 { 765 766 //-------------------------------------------------- 767 // Method04h: M(:,j) is sparse 768 //-------------------------------------------------- 769 770 ASSERT (M_is_sparse || M_is_hyper) ; 771 for ( ; pZ < pZ_end ; pZ++) 772 { 773 int64_t i = Zi [pZ] ; 774 bool mij = false ; 775 int64_t pright = pM_end - 1 ; 776 bool found ; 777 GB_BINARY_SEARCH (i, Mi, pM, pright, found) ; 778 if (found) mij = GB_mcast (Mx, pM, msize) ; 779 if (!mij) GB_COPY_Z ; // mask is complemented 780 } 781 } 782 } 783 784 } 785 else if (zjnz == 0) 786 { 787 788 //---------------------------------------------------------- 789 // Z(:,j) is empty 790 //---------------------------------------------------------- 791 792 if (Mask_comp) 793 { 794 795 //------------------------------------------------------ 796 // mask is complemented 797 //------------------------------------------------------ 798 799 if (mdense) 800 { 801 802 //-------------------------------------------------- 803 // Method04i: M(:,j) is dense 804 //-------------------------------------------------- 805 806 for ( ; pC < pC_end ; pC++) 807 { 808 int64_t i = Ci [pC] ; 809 // mask is dense, lookup M(i,j) 810 pM = pM_first + (i - iM_first) ; 811 ASSERT (i == GBI (Mi, pM, vlen)) ; 812 bool mij = GBB (Mb, pM) && 813 GB_mcast (Mx, pM, msize) ; 814 if (mij) GB_COPY_C ; 815 } 816 817 } 818 else if (cjnz > 32 * mjnz) 819 { 820 821 //-------------------------------------------------- 822 // Method04j: C(:,j) is much denser than M(:,j) 823 //-------------------------------------------------- 824 825 ASSERT (M_is_sparse || M_is_hyper) ; 826 for ( ; pM < pM_end ; pM++) 827 { 828 if (GB_mcast (Mx, pM, msize)) 829 { 830 int64_t i = Mi [pM] ; 831 int64_t pright = pC_end - 1 ; 832 bool found ; 833 GB_BINARY_SEARCH (i, Ci, pC, pright, found); 834 if (found) GB_COPY_C ; 835 } 836 } 837 838 } 839 else if (mjnz > 32 * cjnz) 840 { 841 842 //-------------------------------------------------- 843 // Method04k: M(:,j) is much denser than C(:,j) 844 //-------------------------------------------------- 845 846 ASSERT (M_is_sparse || M_is_hyper) ; 847 for ( ; pC < pC_end ; pC++) 848 { 849 int64_t i = Ci [pC] ; 850 bool mij = false ; 851 int64_t pright = pM_end - 1 ; 852 bool found ; 853 GB_BINARY_SEARCH (i, Mi, pM, pright, found); 854 if (found) mij = GB_mcast (Mx, pM, msize) ; 855 if (mij) GB_COPY_C ; 856 } 857 858 } 859 else 860 { 861 862 //-------------------------------------------------- 863 // Method04l: M(:,j) and C(:,j) about same # entries 864 //-------------------------------------------------- 865 866 ASSERT (M_is_sparse || M_is_hyper) ; 867 while (pM < pM_end && pC < pC_end) 868 { 869 int64_t iM = Mi [pM] ; 870 int64_t i = Ci [pC] ; 871 if (iM < i) 872 { 873 // M(i,j) exists but not C(i,j) 874 pM++ ; 875 } 876 else if (i < iM) 877 { 878 // C(i,j) exists but not M(i,j) 879 pC++ ; 880 } 881 else 882 { 883 // both M(i,j) and C(i,j) exist 884 if (GB_mcast (Mx, pM, msize)) GB_COPY_C ; 885 pM++ ; 886 pC++ ; 887 } 888 } 889 } 890 891 } 892 else 893 { 894 895 //------------------------------------------------------ 896 // non-complemented mask, and Z(:,j) empty 897 //------------------------------------------------------ 898 899 if (mdense) 900 { 901 902 //-------------------------------------------------- 903 // Method04m: M(:,j) is dense 904 //-------------------------------------------------- 905 906 for ( ; pC < pC_end ; pC++) 907 { 908 int64_t i = Ci [pC] ; 909 // mask is dense, lookup M(i,j) 910 pM = pM_first + (i - iM_first) ; 911 ASSERT (i == GBI (Mi, pM, vlen)) ; 912 bool mij = GBB (Mb, pM) && 913 GB_mcast (Mx, pM, msize) ; 914 if (!mij) GB_COPY_C ; 915 } 916 } 917 else 918 { 919 920 //-------------------------------------------------- 921 // Method04n: M(:,j) is sparse 922 //-------------------------------------------------- 923 924 ASSERT (M_is_sparse || M_is_hyper) ; 925 for ( ; pC < pC_end ; pC++) 926 { 927 int64_t i = Ci [pC] ; 928 // M(i,j) false if not present 929 bool mij = false ; 930 int64_t pright = pM_end - 1 ; 931 bool found ; 932 GB_BINARY_SEARCH (i, Mi, pM, pright, found) ; 933 if (found) mij = GB_mcast (Mx, pM, msize) ; 934 if (!mij) GB_COPY_C ; 935 } 936 } 937 } 938 } 939 940 #if defined ( GB_PHASE_2_OF_2 ) 941 ASSERT (pR == pR_end) ; 942 #endif 943 } 944 945 //------------------------------------------------------------------ 946 // final count of nnz (R(:,j)) 947 //------------------------------------------------------------------ 948 949 #if defined ( GB_PHASE_1_OF_2 ) 950 if (fine_task) 951 { 952 TaskList [taskid].pC = rjnz ; 953 } 954 else 955 { 956 Rp [k] = rjnz ; 957 } 958 #endif 959 } 960 } 961 } 962 963