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