1 //------------------------------------------------------------------------------
2 // GB_select: apply a select operator
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 // C<M> = accum (C, select(A,Thunk)) or select(A,Thunk)').
11
12 #define GB_FREE_ALL \
13 { \
14 GB_phbix_free (T) ; \
15 }
16
17 #include "GB_select.h"
18 #include "GB_accum_mask.h"
19
GB_select(GrB_Matrix C,const bool C_replace,const GrB_Matrix M,const bool Mask_comp,const bool Mask_struct,const GrB_BinaryOp accum,const GxB_SelectOp op,const GrB_Matrix A,const GxB_Scalar Thunk_in,const bool A_transpose,GB_Context Context)20 GrB_Info GB_select // C<M> = accum (C, select(A,k)) or select(A',k)
21 (
22 GrB_Matrix C, // input/output matrix for results
23 const bool C_replace, // C descriptor
24 const GrB_Matrix M, // optional mask for C, unused if NULL
25 const bool Mask_comp, // descriptor for M
26 const bool Mask_struct, // if true, use the only structure of M
27 const GrB_BinaryOp accum, // optional accum for Z=accum(C,T)
28 const GxB_SelectOp op, // operator to select the entries
29 const GrB_Matrix A, // input matrix
30 const GxB_Scalar Thunk_in, // optional input for select operator
31 const bool A_transpose, // A matrix descriptor
32 GB_Context Context
33 )
34 {
35
36 //--------------------------------------------------------------------------
37 // check inputs
38 //--------------------------------------------------------------------------
39
40 // C may be aliased with M and/or A
41
42 GB_RETURN_IF_FAULTY_OR_POSITIONAL (accum) ;
43 GB_RETURN_IF_FAULTY (Thunk_in) ;
44 GB_RETURN_IF_NULL_OR_FAULTY (op) ;
45
46 ASSERT_MATRIX_OK (C, "C input for GB_select", GB0) ;
47 ASSERT_MATRIX_OK_OR_NULL (M, "M for GB_select", GB0) ;
48 ASSERT_BINARYOP_OK_OR_NULL (accum, "accum for GB_select", GB0) ;
49 ASSERT_SELECTOP_OK (op, "selectop for GB_select", GB0) ;
50 ASSERT_MATRIX_OK (A, "A input for GB_select", GB0) ;
51 ASSERT_SCALAR_OK_OR_NULL (Thunk_in, "Thunk_in for GB_select", GB0) ;
52
53 struct GB_Matrix_opaque T_header ;
54 GrB_Matrix T = GB_clear_static_header (&T_header) ;
55
56 // check domains and dimensions for C<M> = accum (C,T)
57 GrB_Info info ;
58 GB_OK (GB_compatible (C->type, C, M, Mask_struct, accum, A->type, Context));
59
60 GB_Type_code typecode = A->type->code ;
61 GB_Select_Opcode opcode = op->opcode ;
62
63 // these opcodes are not availabe to the user
64 ASSERT (opcode != GB_RESIZE_opcode) ;
65 ASSERT (opcode != GB_NONZOMBIE_opcode) ;
66
67 // check if the op is a GT, GE, LT, or LE comparator
68 bool op_is_ordered_comparator =
69 opcode == GB_GT_ZERO_opcode || opcode == GB_GT_THUNK_opcode ||
70 opcode == GB_GE_ZERO_opcode || opcode == GB_GE_THUNK_opcode ||
71 opcode == GB_LT_ZERO_opcode || opcode == GB_LT_THUNK_opcode ||
72 opcode == GB_LE_ZERO_opcode || opcode == GB_LE_THUNK_opcode ;
73
74 if (op_is_ordered_comparator)
75 {
76 // built-in GT, GE, LT, and LE operators cannot be used with
77 // user-defined or complex types.
78 if (typecode == GB_UDT_code)
79 {
80 GB_ERROR (GrB_DOMAIN_MISMATCH,
81 "Operator %s not defined for user-defined types", op->name) ;
82 }
83 else if (typecode == GB_FC32_code || typecode == GB_FC64_code)
84 {
85 GB_ERROR (GrB_DOMAIN_MISMATCH,
86 "Operator %s not defined for complex types", op->name) ;
87 }
88 }
89
90 // C = op (A) must be compatible, already checked in GB_compatible
91
92 // A must also be compatible with op->xtype
93 if (!GB_Type_compatible (A->type, op->xtype))
94 {
95 GB_ERROR (GrB_DOMAIN_MISMATCH,
96 "Incompatible type for C=%s(A,Thunk):\n"
97 "input A type [%s]\n"
98 "cannot be typecast to operator input of type [%s]",
99 op->name, A->type->name, op->xtype->name) ;
100 }
101
102 // check the dimensions
103 int64_t tnrows = (A_transpose) ? GB_NCOLS (A) : GB_NROWS (A) ;
104 int64_t tncols = (A_transpose) ? GB_NROWS (A) : GB_NCOLS (A) ;
105 if (GB_NROWS (C) != tnrows || GB_NCOLS (C) != tncols)
106 {
107 GB_ERROR (GrB_DIMENSION_MISMATCH,
108 "Dimensions not compatible:\n"
109 "output is " GBd "-by-" GBd "\n"
110 "input is " GBd "-by-" GBd "%s",
111 GB_NROWS (C), GB_NCOLS (C),
112 tnrows, tncols, A_transpose ? " (transposed)" : "") ;
113 }
114
115 // check if op is (NE, EQ, GT, GE, LT, LE)_THUNK
116 bool op_is_thunk_comparator =
117 (opcode >= GB_NE_THUNK_opcode && opcode <= GB_LE_THUNK_opcode) ;
118
119 // check if op is TRIL, TRIU, DIAG, or OFFDIAG
120 bool op_is_positional = GB_SELECTOP_IS_POSITIONAL (opcode) ;
121
122 // check if op is user-defined
123 bool op_is_user_defined = (opcode >= GB_USER_SELECT_opcode) ;
124
125 int64_t nz_thunk = 0 ;
126 GB_void *restrict xthunk_in = NULL ;
127
128 if (Thunk_in != NULL)
129 {
130
131 // finish any pending work on the Thunk
132 GB_MATRIX_WAIT (Thunk_in) ;
133 nz_thunk = GB_NNZ (Thunk_in) ;
134
135 // if op is TRIL, TRIU, DIAG, or OFFDIAG, Thunk_in must be
136 // compatible with GrB_INT64
137 if (op_is_positional &&
138 !GB_Type_compatible (GrB_INT64, Thunk_in->type))
139 {
140 // Thunk not a built-in type, for a built-in select operator
141 GB_ERROR (GrB_DOMAIN_MISMATCH,
142 "Incompatible type for C=%s(A,Thunk):\n"
143 "input Thunk type [%s]\n"
144 "not compatible with GrB_INT64 input to built-in operator %s",
145 op->name, Thunk_in->type->name, op->name) ;
146 }
147
148 // if op is (NE, EQ, GT, GE, LT, LE)_THUNK, then Thunk must be
149 // compatible with the matrix type
150 if (op_is_thunk_comparator &&
151 !GB_Type_compatible (A->type, Thunk_in->type))
152 {
153 GB_ERROR (GrB_DOMAIN_MISMATCH,
154 "Incompatible type for C=%s(A,Thunk):\n"
155 "input A type [%s] and Thunk type [%s] not compatible",
156 op->name, A->type->name, Thunk_in->type->name) ;
157 }
158
159 // get the pointer to the value of Thunk_in
160 xthunk_in = (GB_void *) Thunk_in->x ;
161 }
162
163 // if op is user-defined, Thunk must match the op->ttype exactly
164 if (op_is_user_defined)
165 {
166 if (op->ttype == NULL && Thunk_in != NULL)
167 {
168 // select operator does not take a Thunk, but one is present
169 GB_ERROR (GrB_DOMAIN_MISMATCH,
170 "User-defined operator %s(A,Thunk) does not take a Thunk\n"
171 "input, but Thunk parameter is non-NULL", op->name) ;
172 }
173 else if (op->ttype != NULL && Thunk_in == NULL)
174 {
175 // select operator takes a Thunk, but Thunk parameter is missing
176 GB_ERROR (GrB_NULL_POINTER,
177 "Required argument is null: [%s]", "Thunk") ;
178 }
179 else if (op->ttype != NULL && Thunk_in != NULL)
180 {
181 // select operator takes a Thunk, and it is present on input.
182 // The types must match exactly.
183 if (op->ttype != Thunk_in->type)
184 {
185 GB_ERROR (GrB_DOMAIN_MISMATCH,
186 "User-defined operator %s(A,Thunk) has a Thunk input\n"
187 "type of [%s], which must exactly match the type of the\n"
188 "Thunk parameter; parameter to GxB_select has type [%s]",
189 op->name, op->ttype->name, Thunk_in->type->name) ;
190 }
191 if (nz_thunk != 1)
192 {
193 GB_ERROR (GrB_INVALID_VALUE,
194 "User-defined operator %s(A,Thunk) has a Thunk input,\n"
195 "which must not be empty", op->name) ;
196 }
197 }
198 }
199
200 // quick return if an empty mask is complemented
201 GB_RETURN_IF_QUICK_MASK (C, C_replace, M, Mask_comp, Mask_struct) ;
202
203 //--------------------------------------------------------------------------
204 // delete any lingering zombies and assemble any pending tuples
205 //--------------------------------------------------------------------------
206
207 GB_MATRIX_WAIT (M) ; // TODO: delay until accum/mask phase
208 GB_MATRIX_WAIT (A) ; // TODO: could tolerate jumbled in some cases
209
210 GB_BURBLE_DENSE (C, "(C %s) ") ;
211 GB_BURBLE_DENSE (M, "(M %s) ") ;
212 GB_BURBLE_DENSE (A, "(A %s) ") ;
213
214 //--------------------------------------------------------------------------
215 // handle the CSR/CSC format and the transposed case
216 //--------------------------------------------------------------------------
217
218 // A and C can be in CSR or CSC format (in any combination), and A can be
219 // transposed first via A_transpose. However, A is not explicitly
220 // transposed first. Instead, the selection operation is modified by
221 // changing the operator, and the resulting matrix T is transposed, if
222 // needed.
223
224 // Instead of explicitly transposing the input matrix A and output T:
225 // If A in CSC format and not transposed: treat as if A and T were CSC
226 // If A in CSC format and transposed: treat as if A and T were CSR
227 // If A in CSR format and not transposed: treat as if A and T were CSR
228 // If A in CSR format and transposed: treat as if A and T were CSC
229
230 bool A_csc = (A->is_csc == !A_transpose) ;
231
232 // The final transpose, if needed, is accomplished in GB_accum_mask, by
233 // tagging T as the same CSR/CSC format as A_csc. If the format of T and C
234 // do not match, GB_accum_mask transposes T, computing C<M>=accum(C,T').
235
236 //--------------------------------------------------------------------------
237 // change the opcode if needed
238 //--------------------------------------------------------------------------
239
240 bool flipij = !A_csc ;
241
242 ASSERT_SCALAR_OK_OR_NULL (Thunk_in, "Thunk_in now GB_select", GB0) ;
243
244 // if A is boolean, get the value of Thunk typecasted to boolean
245 bool bthunk = false ;
246
247 if (typecode == GB_BOOL_code && op_is_thunk_comparator && nz_thunk > 0)
248 {
249 // bthunk = (bool) Thunk_in
250 GB_cast_array ((GB_void *) (&bthunk), GB_BOOL_code, xthunk_in,
251 Thunk_in->type->code, NULL, Thunk_in->type->size, 1, 1) ;
252 }
253
254 int64_t ithunk = 0 ; // ithunk = (int64_t) Thunk (0)
255 bool use_dup = false ;
256 bool is_empty = false ;
257
258 if (op_is_positional)
259 {
260
261 //----------------------------------------------------------------------
262 // tril, triu, diag, offdiag: get k and handle the flip
263 //----------------------------------------------------------------------
264
265 // The built-in operators are modified so they can always work as if A
266 // were in CSC format. If A is not in CSC, then the operation is
267 // flipped.
268 // 0: tril(A,k) becomes triu(A,-k)
269 // 1: triu(A,k) becomes tril(A,-k)
270 // 2: diag(A,k) becomes diag(A,-k)
271 // 3: offdiag(A,k) becomes offdiag(A,-k)
272 // all others Thunk is unchanged
273 // userop(A) row/col indices and dimensions are swapped
274
275 // if Thunk is not present, or has no entries, then k defaults to zero
276 if (nz_thunk > 0)
277 {
278 // ithunk = (int64_t) (Thunk_in (0)) ;
279 GB_cast_array ((GB_void *) &ithunk, GB_INT64_code, xthunk_in,
280 Thunk_in->type->code, NULL, Thunk_in->type->size, 1, 1) ;
281 }
282
283 if (flipij)
284 {
285 ithunk = -ithunk ;
286 if (opcode == GB_TRIL_opcode)
287 {
288 opcode = GB_TRIU_opcode ;
289 }
290 else if (opcode == GB_TRIU_opcode)
291 {
292 opcode = GB_TRIL_opcode ;
293 }
294 flipij = false ;
295 }
296
297 }
298 else
299 {
300
301 //----------------------------------------------------------------------
302 // (NE, EQ, GT, GE, LT, LE) x (0, thunk): handle bool and uint cases
303 //----------------------------------------------------------------------
304
305 switch (opcode)
306 {
307
308 case GB_GT_ZERO_opcode : // A(i,j) > 0
309
310 // bool and uint: rename GxB_GT_ZERO to GxB_NONZERO
311 // user type: return error above
312 switch (typecode)
313 {
314 case GB_BOOL_code :
315 case GB_UINT8_code :
316 case GB_UINT16_code :
317 case GB_UINT32_code :
318 case GB_UINT64_code : opcode = GB_NONZERO_opcode ; break ;
319 default: ;
320 }
321 break ;
322
323 case GB_GE_ZERO_opcode : // A(i,j) >= 0
324
325 // bool and uint: always true; use GB_dup2
326 // user type: return error above
327 switch (typecode)
328 {
329 case GB_BOOL_code :
330 case GB_UINT8_code :
331 case GB_UINT16_code :
332 case GB_UINT32_code :
333 case GB_UINT64_code : use_dup = true ; break ;
334 default: ;
335 }
336 break ;
337
338 case GB_LT_ZERO_opcode : // A(i,j) < 0
339
340 // bool and uint: always false; return an empty matrix
341 // user type: return error above
342 switch (typecode)
343 {
344 case GB_BOOL_code :
345 case GB_UINT8_code :
346 case GB_UINT16_code :
347 case GB_UINT32_code :
348 case GB_UINT64_code : is_empty = true ; break ;
349 default: ;
350 }
351 break ;
352
353 case GB_LE_ZERO_opcode : // A(i,j) <= 0
354
355 // bool and uint: rename GxB_LE_ZERO to GxB_EQ_ZERO
356 // user type: return error above
357 switch (typecode)
358 {
359 case GB_BOOL_code :
360 case GB_UINT8_code :
361 case GB_UINT16_code :
362 case GB_UINT32_code :
363 case GB_UINT64_code : opcode = GB_EQ_ZERO_opcode ; break ;
364 default: ;
365 }
366 break ;
367
368 case GB_NE_THUNK_opcode : // A(i,j) != thunk
369
370 // bool: if thunk is true, rename GxB_NE_THUNK to GxB_EQ_ZERO
371 // if thunk is false, rename GxB_NE_THUNK to GxB_NONZERO
372 if (typecode == GB_BOOL_code)
373 {
374 opcode = (bthunk) ? GB_EQ_ZERO_opcode : GB_NONZERO_opcode ;
375 }
376 break ;
377
378 case GB_EQ_THUNK_opcode : // A(i,j) == thunk
379
380 // bool: if thunk is true, rename GxB_NE_THUNK to GxB_NONZERO
381 // if thunk is false, rename GxB_NE_THUNK to GxB_EQ_ZERO
382 if (typecode == GB_BOOL_code)
383 {
384 opcode = (bthunk) ? GB_NONZERO_opcode : GB_EQ_ZERO_opcode ;
385 }
386 break ;
387
388 case GB_GT_THUNK_opcode : // A(i,j) > thunk
389
390 // bool: if thunk is true, return an empty matrix
391 // if thunk is false, rename GxB_GT_THUNK to GxB_NONZERO
392 // user type: return error above
393 if (typecode == GB_BOOL_code)
394 {
395 if (bthunk)
396 {
397 is_empty = true ;
398 }
399 else
400 {
401 // rename GT_THUNK to NONZERO for boolean
402 opcode = GB_NONZERO_opcode ;
403 }
404 }
405 break ;
406
407 case GB_GE_THUNK_opcode : // A(i,j) >= thunk
408
409 // bool: if thunk is true, rename GxB_GE_THUNK to GxB_NONZERO
410 // if thunk is false, use GB_dup2
411 // user type: return error above
412 if (typecode == GB_BOOL_code)
413 {
414 if (bthunk)
415 {
416 opcode = GB_NONZERO_opcode ;
417 }
418 else
419 {
420 // use dup for GE_THUNK if thunk is false
421 use_dup = true ;
422 }
423 }
424 break ;
425
426 case GB_LT_THUNK_opcode : // A(i,j) < thunk
427
428 // bool: if thunk is true, rename GxB_LT_THUNK to GxB_EQ_ZERO
429 // if thunk is false, return an empty matrix
430 // user type: return error above
431 if (typecode == GB_BOOL_code)
432 {
433 if (bthunk)
434 {
435 opcode = GB_EQ_ZERO_opcode ;
436 }
437 else
438 {
439 // matrix empty for LT_THUNK_BOOL, if thunk false
440 is_empty = true ;
441 }
442 }
443 break ;
444
445 case GB_LE_THUNK_opcode : // A(i,j) <= thunk
446
447 // bool: if thunk is true, use GB_dup2
448 // if thunk is false, rename GxB_LE_ZERO to GxB_EQ_ZERO
449 // user type: return error
450 if (typecode == GB_BOOL_code)
451 {
452 if (bthunk)
453 {
454 // use dup for LE_THUNK if thunk is true
455 use_dup = true ;
456 }
457 else
458 {
459 opcode = GB_EQ_ZERO_opcode ;
460 }
461 }
462 break ;
463
464 default : ; // use the opcode as-is
465 }
466 }
467
468 //--------------------------------------------------------------------------
469 // create T
470 //--------------------------------------------------------------------------
471
472 if (use_dup)
473 {
474 // selectop is always true, so use GB_dup2 to do T = A
475 GB_OK (GB_dup2 (&T, A, true, A->type, Context)) ; // static header
476 }
477 else if (is_empty)
478 {
479 // selectop is always false, so T is an empty matrix
480 GB_OK (GB_new (&T, true, // auto (sparse or hyper), static header
481 A->type, A->vlen, A->vdim, GB_Ap_calloc, A_csc,
482 GxB_SPARSE + GxB_HYPERSPARSE, GB_Global_hyper_switch_get ( ),
483 1, Context)) ;
484 }
485 else
486 {
487 // T = select (A, Thunk)
488 GB_OK (GB_selector (T, opcode, op, flipij, A, ithunk,
489 (op_is_thunk_comparator || op_is_user_defined) ? Thunk_in : NULL,
490 Context)) ;
491 }
492
493 T->is_csc = A_csc ;
494 ASSERT_MATRIX_OK (T, "T=select(A,Thunk) output", GB0) ;
495 ASSERT_MATRIX_OK (C, "C for accum; T=select(A,Thunk) output", GB0) ;
496
497 //--------------------------------------------------------------------------
498 // C<M> = accum (C,T): accumulate the results into C via the mask
499 //--------------------------------------------------------------------------
500
501 return (GB_accum_mask (C, M, NULL, accum, &T, C_replace, Mask_comp,
502 Mask_struct, Context)) ;
503 }
504
505