1 //------------------------------------------------------------------------------
2 // GB_slice_vector: slice a vector for GB_add, GB_emult, and GB_mask
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 // A(:,kA) and B(:,kB) are two long vectors that will be added with GB_add,
11 // GB_emult, or GB_mask, and the work to compute them needs to be split into
12 // multiple tasks. They represent the same vector index j, for:
13
14 // C(:,j) = A(:,j) + B(:,j) in GB_add
15 // C(:,j) = A(:,j) .* B(:,j) in GB_emult
16 // C(:,j)<M(:,j)> = B(:,j) in GB_mask; A is passed in as the input C
17 // union (A->h, B->h) in GB_add_phase0.
18
19 // The vector index j is not needed here. The vectors kA and kB are not
20 // required, either; just the positions where the vectors appear in A and B
21 // (pA_start, pA_end, pB_start, and pB_end).
22
23 // The inputs Mi, Ai, and Bi must be sorted on input.
24
25 // This method finds i so that nnz (A (i:end,kA)) + nnz (B (i:end,kB)) is
26 // roughly equal to target_work. The entries in A(i:end,kA) start at position
27 // pA in Ai and Ax, and the entries in B(i:end,kB) start at position pB in Bi
28 // and Bx. Once the work is split, pM is found for M(i:end,kM), if the mask M
29 // is present.
30
31 // The lists Mi, Ai, and Bi can also be any sorted integer array. This is used
32 // by GB_add_phase0 to construct the set union of A->h and B->h. In this case,
33 // pA_start and pB_start are both zero, and pA_end and pB_end are A->nvec and
34 // B->nvec, respectively.
35
36 // If n = A->vlen = B->vlen, anz = nnz (A (:,kA)), and bnz = nnz (B (:,kB)),
37 // then the total time taken by this function is O(log(n)*(log(anz)+log(bnz))),
38 // or at most O((log(n)^2)).
39
40 // The input matrices M, A, and B are not present here, except for M->i,
41 // A->i, and B->i if they are sparse or hypersparse. They cannot be jumbled.
42 // M, A, and B can have any sparsity structure. If bitmap or full, their
43 // corresponding [A,B,M]->i arrays are NULL.
44
45 #include "GB.h"
46
GB_slice_vector(int64_t * p_i,int64_t * p_pM,int64_t * p_pA,int64_t * p_pB,const int64_t pM_start,const int64_t pM_end,const int64_t * restrict Mi,const int64_t pA_start,const int64_t pA_end,const int64_t * restrict Ai,const int64_t pB_start,const int64_t pB_end,const int64_t * restrict Bi,const int64_t vlen,const double target_work)47 void GB_slice_vector
48 (
49 // output: return i, pA, and pB
50 int64_t *p_i, // work starts at A(i,kA) and B(i,kB)
51 int64_t *p_pM, // M(i:end,kM) starts at pM
52 int64_t *p_pA, // A(i:end,kA) starts at pA
53 int64_t *p_pB, // B(i:end,kB) starts at pB
54 // input:
55 const int64_t pM_start, // M(:,kM) starts at pM_start in Mi,Mx
56 const int64_t pM_end, // M(:,kM) ends at pM_end-1 in Mi,Mx
57 const int64_t *restrict Mi, // indices of M (or NULL)
58 const int64_t pA_start, // A(:,kA) starts at pA_start in Ai,Ax
59 const int64_t pA_end, // A(:,kA) ends at pA_end-1 in Ai,Ax
60 const int64_t *restrict Ai, // indices of A (or NULL)
61 const int64_t pB_start, // B(:,kB) starts at pB_start in Bi,Bx
62 const int64_t pB_end, // B(:,kB) ends at pB_end-1 in Bi,Bx
63 const int64_t *restrict Bi, // indices of B (or NULL)
64 const int64_t vlen, // A->vlen and B->vlen
65 const double target_work // target work
66 )
67 {
68
69 //--------------------------------------------------------------------------
70 // check inputs
71 //--------------------------------------------------------------------------
72
73 ASSERT (p_pA != NULL && p_pB != NULL) ;
74
75 //--------------------------------------------------------------------------
76 // find i, pA, and pB for the start of this task
77 //--------------------------------------------------------------------------
78
79 // search for index i in the range ileft:iright, inclusive
80 int64_t ileft = 0 ;
81 int64_t iright = vlen-1 ;
82 int64_t i = 0 ;
83
84 int64_t aknz = pA_end - pA_start ;
85 int64_t bknz = pB_end - pB_start ;
86 int64_t mknz = pM_end - pM_start ; // zero if M not present
87
88 bool a_empty = (aknz == 0) ;
89 bool b_empty = (bknz == 0) ;
90 bool m_empty = (mknz == 0) ;
91
92 int64_t pM = (m_empty) ? -1 : pM_start ;
93 int64_t pA = (a_empty) ? -1 : pA_start ;
94 int64_t pB = (b_empty) ? -1 : pB_start ;
95
96 while (ileft < iright)
97 {
98
99 //----------------------------------------------------------------------
100 // find the index i in the middle of ileft:iright
101 //----------------------------------------------------------------------
102
103 i = (ileft + iright) >> 1 ;
104
105 //----------------------------------------------------------------------
106 // find where i appears in A(:,kA)
107 //----------------------------------------------------------------------
108
109 if (a_empty)
110 {
111 // Ai is empty so i does not appear
112 pA = -1 ;
113 }
114 else if (aknz == vlen)
115 {
116 // A(:,kA) is dense (bitmap, full, or all entries present)
117 // no need for a binary search
118 pA = pA_start + i ;
119 ASSERT (GBI (Ai, pA, vlen) == i) ;
120 }
121 else
122 {
123 // Ai is an explicit integer list, Ai [pA_start:pA_end-1]
124 ASSERT (aknz > 0) ;
125 pA = pA_start ;
126 bool afound ;
127 int64_t apright = pA_end - 1 ;
128 GB_SPLIT_BINARY_SEARCH (i, Ai, pA, apright, afound) ;
129 ASSERT (GB_IMPLIES (afound, GBI (Ai, pA, vlen) == i)) ;
130 ASSERT (pA_start <= pA && pA <= pA_end) ;
131 }
132
133 ASSERT (GB_IMPLIES (pA > pA_start && pA < pA_end,
134 (GBI (Ai, pA-1, vlen) < i))) ;
135 ASSERT (GB_IMPLIES (pA >= pA_start && pA < pA_end,
136 (GBI (Ai, pA, vlen) >= i ))) ;
137
138 // Ai has been split. If afound is false:
139 // Ai [pA_start : pA-1] < i
140 // Ai [pA : pA_end-1] > i
141 // If afound is true:
142 // Ai [pA_start : pA-1] < i
143 // Ai [pA : pA_end-1] >= i
144 //
145 // in both cases, if i is chosen as the breakpoint, then the
146 // subtask starts at index i, and position pA in Ai,Ax.
147
148 // if A(:,kA) is empty, then pA is -1
149
150 //----------------------------------------------------------------------
151 // find where i appears in B(:,kB)
152 //----------------------------------------------------------------------
153
154 if (b_empty)
155 {
156 // B(:,kB) is empty so i does not appear
157 pB = -1 ;
158 }
159 else if (bknz == vlen)
160 {
161 // B(:,kB) is dense (bitmap, full, or all entries present)
162 // no need for a binary search
163 pB = pB_start + i ;
164 ASSERT (GBI (Bi, pB, vlen) == i) ;
165 }
166 else
167 {
168 // B(:,kB) is sparse, and not empty
169 ASSERT (bknz > 0) ;
170 ASSERT (Bi != NULL) ;
171 pB = pB_start ;
172 bool bfound ;
173 int64_t bpright = pB_end - 1 ;
174 GB_SPLIT_BINARY_SEARCH (i, Bi, pB, bpright, bfound) ;
175 ASSERT (pB_start <= pB && pB <= pB_end) ;
176 }
177 ASSERT (GB_IMPLIES (pB > pB_start && pB < pB_end,
178 (GBI (Bi, pB-1, vlen) < i))) ;
179 ASSERT (GB_IMPLIES (pB >= pB_start && pB < pB_end,
180 (GBI (Bi, pB, vlen) >= i ))) ;
181
182 // Bi has been split. If bfound is false:
183 // Bi [pB_start : pB-1] < i
184 // Bi [pB : pB_end-1] > i
185 // If bfound is true:
186 // Bi [pB_start : pB-1] < i
187 // Bi [pB : pB_end-1] >= i
188 //
189 // in both cases, if i is chosen as the breakpoint, then the
190 // subtask starts at index i, and position pB in Bi,Bx.
191
192 // if B(:,kB) is empty, then pB is -1
193
194 //----------------------------------------------------------------------
195 // determine if the subtask is near the target task size
196 //----------------------------------------------------------------------
197
198 double work = (a_empty ? 0 : (pA_end - pA))
199 + (b_empty ? 0 : (pB_end - pB)) ;
200
201 if (work < 0.9999 * target_work)
202 {
203
204 //------------------------------------------------------------------
205 // work is too low
206 //------------------------------------------------------------------
207
208 // work is too low, so i is too high.
209 // Keep searching in the range (ileft:i), inclusive.
210
211 iright = i ;
212
213 }
214 else if (work > 1.0001 * target_work)
215 {
216
217 //------------------------------------------------------------------
218 // work is too high
219 //------------------------------------------------------------------
220
221 // work is too high, so i is too low.
222 // Keep searching in the range (i+1):iright, inclusive.
223
224 ileft = i + 1 ;
225
226 }
227 else
228 {
229
230 //------------------------------------------------------------------
231 // work is about right; use this result.
232 //------------------------------------------------------------------
233
234 // return i, pA, and pB as the start of this task.
235 ASSERT (0 <= i && i <= vlen) ;
236 ASSERT (pA == -1 || (pA_start <= pA && pA <= pA_end)) ;
237 ASSERT (pB == -1 || (pB_start <= pB && pB <= pB_end)) ;
238 break ;
239 }
240 }
241
242 //--------------------------------------------------------------------------
243 // find where i appears in M(:,kM)
244 //--------------------------------------------------------------------------
245
246 if (m_empty)
247 {
248 pM = -1 ;
249 }
250 else if (mknz == vlen)
251 {
252 // M(:,kM) is dense (bitmap, full, or all entries present)
253 // no need for a binary search
254 pM = pM_start + i ;
255 ASSERT (GBI (Mi, pM, vlen) == i) ;
256 }
257 else
258 {
259 // M(:,kM) is sparse, and not empty
260 ASSERT (mknz > 0) ;
261 ASSERT (Mi != NULL) ;
262 pM = pM_start ;
263 bool mfound ;
264 int64_t mpright = pM_end - 1 ;
265 GB_SPLIT_BINARY_SEARCH (i, Mi, pM, mpright, mfound) ;
266 }
267
268 //--------------------------------------------------------------------------
269 // return result
270 //--------------------------------------------------------------------------
271
272 // pM, pA, and pB partition the three vectors M(:,j), A(:,j), and B(:,j),
273 // or if any vector is empty, their p* pointer is -1.
274
275 ASSERT (GB_IMPLIES ((pM > pM_start && pM < pM_end),
276 GBI (Mi, pM-1, vlen) < i)) ;
277 ASSERT (GB_IMPLIES ((pM >= pM_start && pM < pM_end),
278 GBI (Mi, pM, vlen) >= i)) ;
279 ASSERT (GB_IMPLIES ((pA > pA_start && pA < pA_end),
280 GBI (Ai, pA-1, vlen) < i)) ;
281 ASSERT (GB_IMPLIES ((pA >= pA_start && pA < pA_end),
282 GBI (Ai, pA, vlen) >= i)) ;
283 ASSERT (GB_IMPLIES ((pB > pB_start && pB < pB_end),
284 GBI (Bi, pB-1, vlen) < i)) ;
285 ASSERT (GB_IMPLIES ((pB >= pB_start && pB < pB_end),
286 GBI (Bi, pB, vlen) >= i)) ;
287
288 if (p_i != NULL)
289 {
290 (*p_i) = i ;
291 }
292 if (p_pM != NULL)
293 {
294 (*p_pM) = pM ;
295 }
296 (*p_pA) = pA ;
297 (*p_pB) = pB ;
298 }
299
300