1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
|
//------------------------------------------------------------------------------
// GB_add_sparse_template: C=A+B, C<M>=A+B when C is sparse/hypersparse
//------------------------------------------------------------------------------
// SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2025, All Rights Reserved.
// SPDX-License-Identifier: Apache-2.0
//------------------------------------------------------------------------------
// C is sparse or hypersparse:
// ------------------------------------------
// C = A + B
// ------------------------------------------
// sparse . sparse sparse
// ------------------------------------------
// C <M> = A + B
// ------------------------------------------
// sparse sparse sparse sparse
// sparse sparse sparse bitmap
// sparse sparse sparse full
// sparse sparse bitmap sparse
// sparse sparse bitmap bitmap
// sparse sparse bitmap full
// sparse sparse full sparse
// sparse sparse full bitmap
// sparse sparse full full
// sparse bitmap sparse sparse
// sparse full sparse sparse
// ------------------------------------------
// C <!M> = A + B
// ------------------------------------------
// sparse bitmap sparse sparse
// sparse full sparse sparse
// If all four matrices are sparse/hypersparse, and C<!M>=A+B is being
// computed, then M is passed in as NULL to GB_add_phase*.
// GB_add_sparsity returns apply_mask as false. The methods below do
// not handle the case when C is sparse, M is sparse, and !M is used.
// All other uses of !M when M is sparse result in a bitmap structure
// for C, and this is handled by GB_add_bitmap_template.
// For this case: the mask is done later, so C=A+B is computed here:
// ------------------------------------------
// C <!M> = A + B
// ------------------------------------------
// sparse sparse sparse sparse (mask later)
{
//--------------------------------------------------------------------------
// phase1: count entries in each C(:,j)
// phase2: compute C
//--------------------------------------------------------------------------
#pragma omp parallel for num_threads(C_nthreads) schedule(dynamic,1)
for (taskid = 0 ; taskid < C_ntasks ; taskid++)
{
//----------------------------------------------------------------------
// get the task descriptor
//----------------------------------------------------------------------
int64_t kfirst = TaskList [taskid].kfirst ;
int64_t klast = TaskList [taskid].klast ;
bool fine_task = (klast == -1) ;
int64_t len ;
if (fine_task)
{
// a fine task operates on a slice of a single vector
klast = kfirst ;
len = TaskList [taskid].len ;
}
else
{
// a coarse task operates on one or more whole vectors
len = vlen ;
}
//----------------------------------------------------------------------
// compute all vectors in this task
//----------------------------------------------------------------------
for (int64_t k = kfirst ; k <= klast ; k++)
{
//------------------------------------------------------------------
// get j, the kth vector of C
//------------------------------------------------------------------
int64_t j = GBh_C (Ch, k) ;
#if ( GB_ADD_PHASE == 1 )
int64_t cjnz = 0 ;
#else
int64_t pC, pC_end ;
if (fine_task)
{
// A fine task computes a slice of C(:,j)
pC = TaskList [taskid ].pC ;
pC_end = TaskList [taskid+1].pC ;
ASSERT (GB_IGET (Cp, k) <= pC) ;
ASSERT (pC <= pC_end) ;
ASSERT (pC_end <= GB_IGET (Cp, k+1)) ;
}
else
{
// The vectors of C are never sliced for a coarse task.
pC = GB_IGET (Cp, k ) ;
pC_end = GB_IGET (Cp, k+1) ;
}
int64_t cjnz = pC_end - pC ;
if (cjnz == 0) continue ;
#endif
//------------------------------------------------------------------
// get A(:,j)
//------------------------------------------------------------------
int64_t pA = -1, pA_end = -1 ;
if (fine_task)
{
// A fine task operates on Ai,Ax [pA...pA_end-1], which is
// a subset of the vector A(:,j)
pA = TaskList [taskid].pA ;
pA_end = TaskList [taskid].pA_end ;
}
else
{
// A coarse task operates on the entire vector A (:,j)
int64_t kA = (C_to_A == NULL) ? j : C_to_A [k] ;
if (kA >= 0)
{
pA = GBp_A (Ap, kA, vlen) ;
pA_end = GBp_A (Ap, kA+1, vlen) ;
}
}
int64_t ajnz = pA_end - pA ; // nnz in A(:,j) for this slice
int64_t pA_start = pA ;
bool adense = (ajnz == len) ;
// get the first and last indices in A(:,j) for this vector
int64_t iA_first = -1, iA_last = -1 ;
if (ajnz > 0)
{
iA_first = GBi_A (Ai, pA, vlen) ;
iA_last = GBi_A (Ai, pA_end-1, vlen) ;
}
//------------------------------------------------------------------
// get B(:,j)
//------------------------------------------------------------------
int64_t pB = -1, pB_end = -1 ;
if (fine_task)
{
// A fine task operates on Bi,Bx [pB...pB_end-1], which is
// a subset of the vector B(:,j)
pB = TaskList [taskid].pB ;
pB_end = TaskList [taskid].pB_end ;
}
else
{
// A coarse task operates on the entire vector B (:,j)
int64_t kB = (C_to_B == NULL) ? j : C_to_B [k] ;
if (kB >= 0)
{
pB = GBp_B (Bp, kB, vlen) ;
pB_end = GBp_B (Bp, kB+1, vlen) ;
}
}
int64_t bjnz = pB_end - pB ; // nnz in B(:,j) for this slice
int64_t pB_start = pB ;
bool bdense = (bjnz == len) ;
// get the first and last indices in B(:,j) for this vector
int64_t iB_first = -1, iB_last = -1 ;
if (bjnz > 0)
{
iB_first = GBi_B (Bi, pB, vlen) ;
iB_last = GBi_B (Bi, pB_end-1, vlen) ;
}
//------------------------------------------------------------------
// C(:,j)<optional mask> = A (:,j) + B (:,j) or subvector
//------------------------------------------------------------------
#ifdef GB_JIT_KERNEL
{
#if GB_NO_MASK
{
#include "template/GB_add_sparse_noM.c"
}
#elif (GB_M_IS_SPARSE || GB_M_IS_HYPER)
{
#include "template/GB_add_sparse_M_sparse.c"
}
#else
{
#include "template/GB_add_sparse_M_bitmap.c"
}
#endif
}
#else
{
if (M == NULL)
{
#include "template/GB_add_sparse_noM.c"
}
else if (M_is_sparse_or_hyper)
{
#include "template/GB_add_sparse_M_sparse.c"
}
else
{
#include "template/GB_add_sparse_M_bitmap.c"
}
}
#endif
//------------------------------------------------------------------
// final count of nnz (C (:,j))
//------------------------------------------------------------------
#if ( GB_ADD_PHASE == 1 )
if (fine_task)
{
TaskList [taskid].pC = cjnz ;
}
else
{
GB_ISET (Cp, k, cjnz) ; // Cp [k] = cjnz ;
}
#endif
}
}
}
|
