1 //------------------------------------------------------------------------------
2 // GB_subassign_10_and_18: C(I,J)<M or !M,repl> = A ; using S
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 // Method 10: C(I,J)<M,repl> = A ; using S
11 // Method 18: C(I,J)<!M,repl> = A ; using S
12
13 // M: present
14 // Mask_comp: true or false
15 // C_replace: true
16 // accum: NULL
17 // A: matrix
18 // S: constructed
19
20 // C: not bitmap: use GB_bitmap_assign instead
21 // M, A: any sparsity structure.
22
23 #include "GB_subassign_methods.h"
24
GB_subassign_10_and_18(GrB_Matrix C,const GrB_Index * I,const int64_t ni,const int64_t nI,const int Ikind,const int64_t Icolon[3],const GrB_Index * J,const int64_t nj,const int64_t nJ,const int Jkind,const int64_t Jcolon[3],const GrB_Matrix M,const bool Mask_struct,const bool Mask_comp,const GrB_Matrix A,GB_Context Context)25 GrB_Info GB_subassign_10_and_18
26 (
27 GrB_Matrix C,
28 // input:
29 const GrB_Index *I,
30 const int64_t ni,
31 const int64_t nI,
32 const int Ikind,
33 const int64_t Icolon [3],
34 const GrB_Index *J,
35 const int64_t nj,
36 const int64_t nJ,
37 const int Jkind,
38 const int64_t Jcolon [3],
39 const GrB_Matrix M,
40 const bool Mask_struct, // if true, use the only structure of M
41 const bool Mask_comp, // if true, !M, else use M
42 const GrB_Matrix A,
43 GB_Context Context
44 )
45 {
46
47 //--------------------------------------------------------------------------
48 // check inputs
49 //--------------------------------------------------------------------------
50
51 ASSERT (!GB_IS_BITMAP (C)) ; ASSERT (!GB_IS_FULL (C)) ;
52 ASSERT (!GB_aliased (C, M)) ; // NO ALIAS of C==M
53 ASSERT (!GB_aliased (C, A)) ; // NO ALIAS of C==A
54
55 //--------------------------------------------------------------------------
56 // S = C(I,J)
57 //--------------------------------------------------------------------------
58
59 GB_EMPTY_TASKLIST ;
60 GB_OK (GB_subassign_symbolic (S, C, I, ni, J, nj, true, Context)) ;
61
62 //--------------------------------------------------------------------------
63 // get inputs
64 //--------------------------------------------------------------------------
65
66 GB_MATRIX_WAIT_IF_JUMBLED (M) ;
67 GB_MATRIX_WAIT_IF_JUMBLED (A) ;
68
69 GB_GET_C ; // C must not be bitmap
70 GB_GET_MASK ;
71 GB_GET_A ;
72 GB_GET_S ;
73 GrB_BinaryOp accum = NULL ;
74
75 //--------------------------------------------------------------------------
76 // Method 10: C(I,J)<M,repl> = A ; using S
77 // Method 18: C(I,J)<!M,repl> = A ; using S
78 //--------------------------------------------------------------------------
79
80 // Time: Optimal. Omega (nnz(A)+nnz(S)), since all entries in S+A must be
81 // traversed, and the corresponding entry in M (even if not present)
82 // determines the action to take. M can add a log(m) factor if sparse.
83
84 //--------------------------------------------------------------------------
85 // Parallel: A+S (Methods 02, 04, 09, 10, 11, 12, 14, 16, 18, 20)
86 //--------------------------------------------------------------------------
87
88 if (A_is_bitmap)
89 {
90 // all of IxJ must be examined
91 GB_SUBASSIGN_IXJ_SLICE ;
92 }
93 else
94 {
95 // traverse all A+S
96 GB_SUBASSIGN_TWO_SLICE (A, S) ;
97 }
98
99 //--------------------------------------------------------------------------
100 // phase 1: create zombies, update entries, and count pending tuples
101 //--------------------------------------------------------------------------
102
103 if (A_is_bitmap)
104 {
105
106 //----------------------------------------------------------------------
107 // phase1: A is bitmap
108 //----------------------------------------------------------------------
109
110 #pragma omp parallel for num_threads(nthreads) schedule(dynamic,1) \
111 reduction(+:nzombies)
112 for (taskid = 0 ; taskid < ntasks ; taskid++)
113 {
114
115 //------------------------------------------------------------------
116 // get the task descriptor
117 //------------------------------------------------------------------
118
119 GB_GET_IXJ_TASK_DESCRIPTOR_PHASE1 (iA_start, iA_end) ;
120
121 //------------------------------------------------------------------
122 // compute all vectors in this task
123 //------------------------------------------------------------------
124
125 for (int64_t j = kfirst ; j <= klast ; j++)
126 {
127
128 //--------------------------------------------------------------
129 // get S(iA_start:iA_end,j)
130 //--------------------------------------------------------------
131
132 GB_GET_VECTOR_FOR_IXJ (S, iA_start) ;
133 int64_t pA_start = j * Avlen ;
134
135 //--------------------------------------------------------------
136 // get M(:,j)
137 //--------------------------------------------------------------
138
139 int64_t pM_start, pM_end ;
140 GB_VECTOR_LOOKUP (pM_start, pM_end, M, j) ;
141 bool mjdense = (pM_end - pM_start) == Mvlen ;
142
143 //--------------------------------------------------------------
144 // do a 2-way merge of S(iA_start:iA_end,j) and A(ditto,j)
145 //--------------------------------------------------------------
146
147 for (int64_t iA = iA_start ; iA < iA_end ; iA++)
148 {
149 int64_t pA = pA_start + iA ;
150 bool Sfound = (pS < pS_end) && (GBI (Si, pS, Svlen) == iA) ;
151 bool Afound = Ab [pA] ;
152
153 if (Sfound && !Afound)
154 {
155 // S (i,j) is present but A (i,j) is not
156 // ----[C . 1] or [X . 1]-------------------------------
157 // [C . 1]: action: ( delete ): becomes zombie
158 // [X . 1]: action: ( X ): still zombie
159 // ----[C . 0] or [X . 0]-------------------------------
160 // [X . 0]: action: ( X ): still a zombie
161 // [C . 0]: C_repl: action: ( delete ): becomes zombie
162 GB_C_S_LOOKUP ;
163 GB_DELETE_ENTRY ;
164 GB_NEXT (S) ;
165 }
166 else if (!Sfound && Afound)
167 {
168 // S (i,j) is not present, A (i,j) is present
169 GB_MIJ_BINARY_SEARCH_OR_DENSE_LOOKUP (iA) ;
170 if (Mask_comp) mij = !mij ;
171 if (mij)
172 {
173 // ----[. A 1]--------------------------------------
174 // [. A 1]: action: ( insert )
175 task_pending++ ;
176 }
177 }
178 else if (Sfound && Afound)
179 {
180 // both S (i,j) and A (i,j) present
181 GB_MIJ_BINARY_SEARCH_OR_DENSE_LOOKUP (iA) ;
182 if (Mask_comp) mij = !mij ;
183 GB_C_S_LOOKUP ;
184 if (mij)
185 {
186 // ----[C A 1] or [X A 1]---------------------------
187 // [C A 1]: action: ( =A ): A to C no accum
188 // [X A 1]: action: ( undelete ): zombie lives
189 GB_noaccum_C_A_1_matrix ;
190 }
191 else
192 {
193 // ----[C A 0] or [X A 0]---------------------------
194 // [X A 0]: action: ( X ): still a zombie
195 // [C A 0]: C_repl: action: ( delete ): now zombie
196 GB_DELETE_ENTRY ;
197 }
198 GB_NEXT (S) ;
199 }
200 }
201 }
202 GB_PHASE1_TASK_WRAPUP ;
203 }
204
205 }
206 else
207 {
208
209 //----------------------------------------------------------------------
210 // phase1: A is hypersparse, sparse, or full
211 //----------------------------------------------------------------------
212
213 #pragma omp parallel for num_threads(nthreads) schedule(dynamic,1) \
214 reduction(+:nzombies)
215 for (taskid = 0 ; taskid < ntasks ; taskid++)
216 {
217
218 //------------------------------------------------------------------
219 // get the task descriptor
220 //------------------------------------------------------------------
221
222 GB_GET_TASK_DESCRIPTOR_PHASE1 ;
223
224 //------------------------------------------------------------------
225 // compute all vectors in this task
226 //------------------------------------------------------------------
227
228 for (int64_t k = kfirst ; k <= klast ; k++)
229 {
230
231 //--------------------------------------------------------------
232 // get A(:,j) and S(:,j)
233 //--------------------------------------------------------------
234
235 int64_t j = GBH (Zh, k) ;
236 GB_GET_MAPPED (pA, pA_end, pA, pA_end, Ap, j, k, Z_to_X, Avlen);
237 GB_GET_MAPPED (pS, pS_end, pB, pB_end, Sp, j, k, Z_to_S, Svlen);
238
239 //--------------------------------------------------------------
240 // get M(:,j)
241 //--------------------------------------------------------------
242
243 int64_t pM_start, pM_end ;
244 GB_VECTOR_LOOKUP (pM_start, pM_end, M, j) ;
245 bool mjdense = (pM_end - pM_start) == Mvlen ;
246
247 //--------------------------------------------------------------
248 // do a 2-way merge of S(:,j) and A(:,j)
249 //--------------------------------------------------------------
250
251 // jC = J [j] ; or J is a colon expression
252 // int64_t jC = GB_ijlist (J, j, Jkind, Jcolon) ;
253
254 // while both list S (:,j) and A (:,j) have entries
255 while (pS < pS_end && pA < pA_end)
256 {
257 int64_t iS = GBI (Si, pS, Svlen) ;
258 int64_t iA = GBI (Ai, pA, Avlen) ;
259
260 if (iS < iA)
261 {
262 // S (i,j) is present but A (i,j) is not
263 // ----[C . 1] or [X . 1]-------------------------------
264 // [C . 1]: action: ( delete ): becomes zombie
265 // [X . 1]: action: ( X ): still zombie
266 // ----[C . 0] or [X . 0]-------------------------------
267 // [X . 0]: action: ( X ): still a zombie
268 // [C . 0]: C_repl: action: ( delete ): becomes zombie
269 GB_C_S_LOOKUP ;
270 GB_DELETE_ENTRY ;
271 GB_NEXT (S) ;
272 }
273 else if (iA < iS)
274 {
275 // S (i,j) is not present, A (i,j) is present
276 GB_MIJ_BINARY_SEARCH_OR_DENSE_LOOKUP (iA) ;
277 if (Mask_comp) mij = !mij ;
278 if (mij)
279 {
280 // ----[. A 1]--------------------------------------
281 // [. A 1]: action: ( insert )
282 task_pending++ ;
283 }
284 GB_NEXT (A) ;
285 }
286 else
287 {
288 // both S (i,j) and A (i,j) present
289 GB_MIJ_BINARY_SEARCH_OR_DENSE_LOOKUP (iA) ;
290 if (Mask_comp) mij = !mij ;
291 GB_C_S_LOOKUP ;
292 if (mij)
293 {
294 // ----[C A 1] or [X A 1]---------------------------
295 // [C A 1]: action: ( =A ): A to C no accum
296 // [X A 1]: action: ( undelete ): zombie lives
297 GB_noaccum_C_A_1_matrix ;
298 }
299 else
300 {
301 // ----[C A 0] or [X A 0]---------------------------
302 // [X A 0]: action: ( X ): still a zombie
303 // [C A 0]: C_repl: action: ( delete ): now zombie
304 GB_DELETE_ENTRY ;
305 }
306 GB_NEXT (S) ;
307 GB_NEXT (A) ;
308 }
309 }
310
311 // while list S (:,j) has entries. List A (:,j) exhausted.
312 while (pS < pS_end)
313 {
314 // ----[C . 1] or [X . 1]-----------------------------------
315 // S (i,j) is present but A (i,j) is not
316 // [C . 1]: action: ( delete ): becomes zombie
317 // [X . 1]: action: ( X ): still a zombie
318 GB_C_S_LOOKUP ;
319 GB_DELETE_ENTRY ;
320 GB_NEXT (S) ;
321 }
322
323 // while list A (:,j) has entries. List S (:,j) exhausted.
324 while (pA < pA_end)
325 {
326 // S (i,j) is not present, A (i,j) is present
327 int64_t iA = GBI (Ai, pA, Avlen) ;
328 GB_MIJ_BINARY_SEARCH_OR_DENSE_LOOKUP (iA) ;
329 if (Mask_comp) mij = !mij ;
330 if (mij)
331 {
332 // ----[. A 1]------------------------------------------
333 // [. A 1]: action: ( insert )
334 task_pending++ ;
335 }
336 GB_NEXT (A) ;
337 }
338 }
339
340 GB_PHASE1_TASK_WRAPUP ;
341 }
342 }
343
344 //--------------------------------------------------------------------------
345 // phase 2: insert pending tuples
346 //--------------------------------------------------------------------------
347
348 GB_PENDING_CUMSUM ;
349
350 if (A_is_bitmap)
351 {
352
353 //----------------------------------------------------------------------
354 // phase2: A is bitmap
355 //----------------------------------------------------------------------
356
357 #pragma omp parallel for num_threads(nthreads) schedule(dynamic,1) \
358 reduction(&&:pending_sorted)
359 for (taskid = 0 ; taskid < ntasks ; taskid++)
360 {
361
362 //------------------------------------------------------------------
363 // get the task descriptor
364 //------------------------------------------------------------------
365
366 GB_GET_IXJ_TASK_DESCRIPTOR_PHASE2 (iA_start, iA_end) ;
367
368 //------------------------------------------------------------------
369 // compute all vectors in this task
370 //------------------------------------------------------------------
371
372 for (int64_t j = kfirst ; j <= klast ; j++)
373 {
374
375 //--------------------------------------------------------------
376 // get S(iA_start:iA_end,j)
377 //--------------------------------------------------------------
378
379 GB_GET_VECTOR_FOR_IXJ (S, iA_start) ;
380 int64_t pA_start = j * Avlen ;
381
382 //--------------------------------------------------------------
383 // get M(:,j)
384 //--------------------------------------------------------------
385
386 int64_t pM_start, pM_end ;
387 GB_VECTOR_LOOKUP (pM_start, pM_end, M, j) ;
388 bool mjdense = (pM_end - pM_start) == Mvlen ;
389
390 //--------------------------------------------------------------
391 // do a 2-way merge of S(iA_start:iA_end,j) and A(ditto,j)
392 //--------------------------------------------------------------
393
394 // jC = J [j] ; or J is a colon expression
395 int64_t jC = GB_ijlist (J, j, Jkind, Jcolon) ;
396
397 for (int64_t iA = iA_start ; iA < iA_end ; iA++)
398 {
399 int64_t pA = pA_start + iA ;
400 bool Sfound = (pS < pS_end) && (GBI (Si, pS, Svlen) == iA) ;
401 bool Afound = Ab [pA] ;
402 if (!Sfound && Afound)
403 {
404 // S (i,j) is not present, A (i,j) is present
405 GB_MIJ_BINARY_SEARCH_OR_DENSE_LOOKUP (iA) ;
406 if (Mask_comp) mij = !mij ;
407 if (mij)
408 {
409 // ----[. A 1]--------------------------------------
410 // [. A 1]: action: ( insert )
411 int64_t iC = GB_ijlist (I, iA, Ikind, Icolon) ;
412 GB_PENDING_INSERT (Ax +(pA*asize)) ;
413 }
414 }
415 else if (Sfound)
416 {
417 // S (i,j) present
418 GB_NEXT (S) ;
419 }
420 }
421 }
422 GB_PHASE2_TASK_WRAPUP ;
423 }
424
425 }
426 else
427 {
428
429 //----------------------------------------------------------------------
430 // phase2: A is hypersparse, sparse, or full
431 //----------------------------------------------------------------------
432
433 #pragma omp parallel for num_threads(nthreads) schedule(dynamic,1) \
434 reduction(&&:pending_sorted)
435 for (taskid = 0 ; taskid < ntasks ; taskid++)
436 {
437
438 //------------------------------------------------------------------
439 // get the task descriptor
440 //------------------------------------------------------------------
441
442 GB_GET_TASK_DESCRIPTOR_PHASE2 ;
443
444 //------------------------------------------------------------------
445 // compute all vectors in this task
446 //------------------------------------------------------------------
447
448 for (int64_t k = kfirst ; k <= klast ; k++)
449 {
450
451 //--------------------------------------------------------------
452 // get A(:,j) and S(:,j)
453 //--------------------------------------------------------------
454
455 int64_t j = GBH (Zh, k) ;
456 GB_GET_MAPPED (pA, pA_end, pA, pA_end, Ap, j, k, Z_to_X, Avlen);
457 GB_GET_MAPPED (pS, pS_end, pB, pB_end, Sp, j, k, Z_to_S, Svlen);
458
459 //--------------------------------------------------------------
460 // get M(:,j)
461 //--------------------------------------------------------------
462
463 int64_t pM_start, pM_end ;
464 GB_VECTOR_LOOKUP (pM_start, pM_end, M, j) ;
465 bool mjdense = (pM_end - pM_start) == Mvlen ;
466
467 //--------------------------------------------------------------
468 // do a 2-way merge of S(:,j) and A(:,j)
469 //--------------------------------------------------------------
470
471 // jC = J [j] ; or J is a colon expression
472 int64_t jC = GB_ijlist (J, j, Jkind, Jcolon) ;
473
474 // while both list S (:,j) and A (:,j) have entries
475 while (pS < pS_end && pA < pA_end)
476 {
477 int64_t iS = GBI (Si, pS, Svlen) ;
478 int64_t iA = GBI (Ai, pA, Avlen) ;
479
480 if (iS < iA)
481 {
482 // S (i,j) is present but A (i,j) is not
483 GB_NEXT (S) ;
484 }
485 else if (iA < iS)
486 {
487 // S (i,j) is not present, A (i,j) is present
488 GB_MIJ_BINARY_SEARCH_OR_DENSE_LOOKUP (iA) ;
489 if (Mask_comp) mij = !mij ;
490 if (mij)
491 {
492 // ----[. A 1]--------------------------------------
493 // [. A 1]: action: ( insert )
494 int64_t iC = GB_ijlist (I, iA, Ikind, Icolon) ;
495 GB_PENDING_INSERT (Ax +(pA*asize)) ;
496 }
497 GB_NEXT (A) ;
498 }
499 else
500 {
501 // both S (i,j) and A (i,j) present
502 GB_NEXT (S) ;
503 GB_NEXT (A) ;
504 }
505 }
506
507 // while list A (:,j) has entries. List S (:,j) exhausted.
508 while (pA < pA_end)
509 {
510 // S (i,j) is not present, A (i,j) is present
511 int64_t iA = GBI (Ai, pA, Avlen) ;
512 GB_MIJ_BINARY_SEARCH_OR_DENSE_LOOKUP (iA) ;
513 if (Mask_comp) mij = !mij ;
514 if (mij)
515 {
516 // ----[. A 1]------------------------------------------
517 // [. A 1]: action: ( insert )
518 int64_t iC = GB_ijlist (I, iA, Ikind, Icolon) ;
519 GB_PENDING_INSERT (Ax +(pA*asize)) ;
520 }
521 GB_NEXT (A) ;
522 }
523 }
524
525 GB_PHASE2_TASK_WRAPUP ;
526 }
527 }
528
529 //--------------------------------------------------------------------------
530 // finalize the matrix and return result
531 //--------------------------------------------------------------------------
532
533 GB_SUBASSIGN_WRAPUP ;
534 }
535
536