1 //------------------------------------------------------------------------------
2 // gbselect: select entries from a GraphBLAS matrix
3 //------------------------------------------------------------------------------
4
5 // SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2021, All Rights Reserved.
6 // SPDX-License-Identifier: GPL-3.0-or-later
7
8 //------------------------------------------------------------------------------
9
10 // gbselect is an interface to GxB_Matrix_select.
11
12 // Usage:
13
14 // C = gbselect (op, A)
15 // C = gbselect (op, A, desc)
16 // C = gbselect (op, A, b, desc)
17
18 // C = gbselect (Cin, accum, op, A, desc)
19 // C = gbselect (Cin, accum, op, A, b, desc)
20
21 // C = gbselect (Cin, M, op, A, desc)
22 // C = gbselect (Cin, M, op, A, b, desc)
23
24 // C = gbselect (Cin, M, accum, op, A, desc)
25 // C = gbselect (Cin, M, accum, op, A, b, desc)
26
27 // If Cin is not present then it is implicitly a matrix with no entries, of the
28 // right size (which depends on A, and the descriptor). The type if Cin, if
29 // not present, is determined by the ztype of the accum, if present, or
30 // otherwise it has the same time as A.
31
32 // If op is '==' or '~=' and b is a NaN, and A has type GrB_FP32, GrB_FP64,
33 // GxB_FC32, or GxB_FC64, then a user-defined operator is used instead of
34 // GxB_EQ_THUNK or GxB_NE_THUNK.
35
36 // The 'tril', 'triu', 'diag', 'offdiag', and 2-input operators all require
37 // the b scalar. The b scalar must not appear for the '*0' operators.
38
39 #include "gb_matlab.h"
40
41 #define USAGE "usage: C = GrB.select (Cin, M, accum, op, A, b, desc)"
42
43 //------------------------------------------------------------------------------
44 // nan operators
45 //------------------------------------------------------------------------------
46
gb_isnan32(GrB_Index i,GrB_Index j,const void * x,const void * b)47 bool gb_isnan32 (GrB_Index i, GrB_Index j, const void *x, const void *b)
48 {
49 float aij = * ((float *) x) ;
50 return (isnan (aij)) ;
51 }
52
gb_isnan64(GrB_Index i,GrB_Index j,const void * x,const void * b)53 bool gb_isnan64 (GrB_Index i, GrB_Index j, const void *x, const void *b)
54 {
55 double aij = * ((double *) x) ;
56 return (isnan (aij)) ;
57 }
58
gb_isnotnan32(GrB_Index i,GrB_Index j,const void * x,const void * b)59 bool gb_isnotnan32 (GrB_Index i, GrB_Index j, const void *x, const void *b)
60 {
61 float aij = * ((float *) x) ;
62 return (!isnan (aij)) ;
63 }
64
gb_isnotnan64(GrB_Index i,GrB_Index j,const void * x,const void * b)65 bool gb_isnotnan64 (GrB_Index i, GrB_Index j, const void *x, const void *b)
66 {
67 double aij = * ((double *) x) ;
68 return (!isnan (aij)) ;
69 }
70
gb_isnanfc32(GrB_Index i,GrB_Index j,const void * x,const void * b)71 bool gb_isnanfc32 (GrB_Index i, GrB_Index j, const void *x, const void *b)
72 {
73 GxB_FC32_t aij = * ((GxB_FC32_t *) x) ;
74 return (isnan (crealf (aij)) || isnan (cimagf (aij))) ;
75 }
76
gb_isnanfc64(GrB_Index i,GrB_Index j,const void * x,const void * b)77 bool gb_isnanfc64 (GrB_Index i, GrB_Index j, const void *x, const void *b)
78 {
79 GxB_FC64_t aij = * ((GxB_FC64_t *) x) ;
80 return (isnan (creal (aij)) || isnan (cimag (aij))) ;
81 }
82
gb_isnotnanfc32(GrB_Index i,GrB_Index j,const void * x,const void * b)83 bool gb_isnotnanfc32 (GrB_Index i, GrB_Index j, const void *x, const void *b)
84 {
85 GxB_FC32_t aij = * ((GxB_FC32_t *) x) ;
86 return (!isnan (crealf (aij)) && !isnan (cimagf (aij))) ;
87 }
88
gb_isnotnanfc64(GrB_Index i,GrB_Index j,const void * x,const void * b)89 bool gb_isnotnanfc64 (GrB_Index i, GrB_Index j, const void *x, const void *b)
90 {
91 GxB_FC64_t aij = * ((GxB_FC64_t *) x) ;
92 return (!isnan (creal (aij)) && !isnan (cimag (aij))) ;
93 }
94
95 //------------------------------------------------------------------------------
96 // gbselect mexFunction
97 //------------------------------------------------------------------------------
98
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])99 void mexFunction
100 (
101 int nargout,
102 mxArray *pargout [ ],
103 int nargin,
104 const mxArray *pargin [ ]
105 )
106 {
107
108 //--------------------------------------------------------------------------
109 // check inputs
110 //--------------------------------------------------------------------------
111
112 gb_usage (nargin >= 2 && nargin <= 7 && nargout <= 2, USAGE) ;
113
114 //--------------------------------------------------------------------------
115 // find the arguments
116 //--------------------------------------------------------------------------
117
118 mxArray *Matrix [4], *String [2], *Cell [2] ;
119 base_enum_t base ;
120 kind_enum_t kind ;
121 GxB_Format_Value fmt ;
122 int nmatrices, nstrings, ncells, sparsity ;
123 GrB_Descriptor desc ;
124 gb_get_mxargs (nargin, pargin, USAGE, Matrix, &nmatrices, String, &nstrings,
125 Cell, &ncells, &desc, &base, &kind, &fmt, &sparsity) ;
126
127 CHECK_ERROR (nmatrices < 1 || nstrings < 1 || ncells > 0, USAGE) ;
128
129 //--------------------------------------------------------------------------
130 // get the select operator
131 //--------------------------------------------------------------------------
132
133 GxB_SelectOp op = gb_mxstring_to_selectop (String [nstrings-1]) ;
134 bool b_required =
135 (op == GxB_TRIL) || (op == GxB_TRIU) ||
136 (op == GxB_DIAG) || (op == GxB_OFFDIAG) ||
137 (op == GxB_NE_THUNK) || (op == GxB_EQ_THUNK) ||
138 (op == GxB_GT_THUNK) || (op == GxB_GE_THUNK) ||
139 (op == GxB_LT_THUNK) || (op == GxB_LE_THUNK) ;
140
141 //--------------------------------------------------------------------------
142 // get the matrices
143 //--------------------------------------------------------------------------
144
145 GrB_Type atype, ctype = NULL ;
146 GrB_Matrix C = NULL, M = NULL, A, b = NULL ;
147
148 if (b_required)
149 {
150 if (nmatrices == 1)
151 {
152 ERROR ("select operator input is missing") ;
153 }
154 else if (nmatrices == 2)
155 {
156 A = gb_get_shallow (Matrix [0]) ;
157 b = gb_get_shallow (Matrix [1]) ;
158 }
159 else if (nmatrices == 3)
160 {
161 C = gb_get_deep (Matrix [0]) ;
162 A = gb_get_shallow (Matrix [1]) ;
163 b = gb_get_shallow (Matrix [2]) ;
164 }
165 else // if (nmatrices == 4)
166 {
167 C = gb_get_deep (Matrix [0]) ;
168 M = gb_get_shallow (Matrix [1]) ;
169 A = gb_get_shallow (Matrix [2]) ;
170 b = gb_get_shallow (Matrix [3]) ;
171 }
172 }
173 else
174 {
175 if (nmatrices == 1)
176 {
177 A = gb_get_shallow (Matrix [0]) ;
178 }
179 else if (nmatrices == 2)
180 {
181 C = gb_get_deep (Matrix [0]) ;
182 A = gb_get_shallow (Matrix [1]) ;
183 }
184 else if (nmatrices == 3)
185 {
186 C = gb_get_deep (Matrix [0]) ;
187 M = gb_get_shallow (Matrix [1]) ;
188 A = gb_get_shallow (Matrix [2]) ;
189 }
190 else // if (nmatrices == 4)
191 {
192 ERROR (USAGE) ;
193 }
194 }
195
196 OK (GxB_Matrix_type (&atype, A)) ;
197 if (C != NULL)
198 {
199 OK (GxB_Matrix_type (&ctype, C)) ;
200 }
201
202 //--------------------------------------------------------------------------
203 // get the accum operator
204 //--------------------------------------------------------------------------
205
206 GrB_BinaryOp accum = NULL ;
207
208 if (nstrings > 1)
209 {
210 // if accum appears, then Cin must also appear
211 CHECK_ERROR (C == NULL, USAGE) ;
212 accum = gb_mxstring_to_binop (String [0], ctype, ctype) ;
213 }
214
215 //--------------------------------------------------------------------------
216 // construct C if not present on input
217 //--------------------------------------------------------------------------
218
219 // If C is NULL, then it is not present on input.
220 // Construct C of the right size and type.
221
222 if (C == NULL)
223 {
224 // get the descriptor contents to determine if A is transposed
225 GrB_Desc_Value in0 ;
226 OK (GxB_Desc_get (desc, GrB_INP0, &in0)) ;
227 bool A_transpose = (in0 == GrB_TRAN) ;
228
229 // get the size of A
230 GrB_Index anrows, ancols ;
231 OK (GrB_Matrix_nrows (&anrows, A)) ;
232 OK (GrB_Matrix_ncols (&ancols, A)) ;
233
234 // determine the size of C
235 GrB_Index cnrows = (A_transpose) ? ancols : anrows ;
236 GrB_Index cncols = (A_transpose) ? anrows : ancols ;
237
238 // C has the same type as A
239 OK (GxB_Matrix_type (&ctype, A)) ;
240
241 // create the matrix C and set its format and sparsity
242 fmt = gb_get_format (cnrows, cncols, A, NULL, fmt) ;
243 sparsity = gb_get_sparsity (A, NULL, sparsity) ;
244 C = gb_new (ctype, cnrows, cncols, fmt, sparsity) ;
245 }
246
247 //--------------------------------------------------------------------------
248 // handle the NaN case
249 //--------------------------------------------------------------------------
250
251 GrB_BinaryOp nan_test = NULL ;
252 GrB_Matrix b2 = b ;
253
254 if (b != NULL)
255 {
256 // check if b is NaN
257 GrB_Type b_type ;
258 OK (GxB_Matrix_type (&b_type, b)) ;
259 bool b_is_nan = false ;
260 if (b_type == GrB_FP32)
261 {
262 float b_value = 0 ;
263 OK (GrB_Matrix_extractElement_FP32 (&b_value, b, 0, 0)) ;
264 b_is_nan = isnan (b_value) ;
265 }
266 else if (b_type == GrB_FP64)
267 {
268 double b_value = 0 ;
269 OK (GrB_Matrix_extractElement_FP64 (&b_value, b, 0, 0)) ;
270 b_is_nan = isnan (b_value) ;
271 }
272 else if (b_type == GxB_FC32)
273 {
274 GxB_FC32_t b_value = GxB_CMPLXF (0, 0) ;
275 OK (GxB_Matrix_extractElement_FC32 (&b_value, b, 0, 0)) ;
276 b_is_nan = GB_cisnanf (b_value) ;
277 }
278 else if (b_type == GxB_FC64)
279 {
280 GxB_FC64_t b_value = GxB_CMPLX (0, 0) ;
281 OK (GxB_Matrix_extractElement_FC64 (&b_value, b, 0, 0)) ;
282 b_is_nan = GB_cisnan (b_value) ;
283 }
284
285 if (b_is_nan)
286 {
287 // b is NaN; create a new nan_test operator if it should be used
288 // instead of the built-in GxB_EQ_THUNK or GxB_NE_THUNK operators.
289 // These operators do not need a b input, since it is now known
290 // to be a NaN.
291
292 if (op == GxB_EQ_THUNK)
293 {
294 if (atype == GrB_FP32)
295 {
296 OK (GxB_SelectOp_new (&nan_test, gb_isnan32,
297 GrB_FP32, NULL)) ;
298 }
299 else if (atype == GrB_FP64)
300 {
301 OK (GxB_SelectOp_new (&nan_test, gb_isnan64,
302 GrB_FP64, NULL)) ;
303 }
304 else if (atype == GxB_FC32)
305 {
306 OK (GxB_SelectOp_new (&nan_test, gb_isnanfc32,
307 GxB_FC32, NULL)) ;
308 }
309 else if (atype == GxB_FC64)
310 {
311 OK (GxB_SelectOp_new (&nan_test, gb_isnanfc64,
312 GxB_FC64, NULL)) ;
313 }
314 }
315 else if (op == GxB_NE_THUNK)
316 {
317 if (atype == GrB_FP32)
318 {
319 OK (GxB_SelectOp_new (&nan_test, gb_isnotnan32,
320 GrB_FP32, NULL)) ;
321 }
322 else if (atype == GrB_FP64)
323 {
324 OK (GxB_SelectOp_new (&nan_test, gb_isnotnan64,
325 GrB_FP64, NULL)) ;
326 }
327 else if (atype == GxB_FC32)
328 {
329 OK (GxB_SelectOp_new (&nan_test, gb_isnotnanfc32,
330 GxB_FC32, NULL)) ;
331 }
332 else if (atype == GxB_FC64)
333 {
334 OK (GxB_SelectOp_new (&nan_test, gb_isnotnanfc64,
335 GxB_FC64, NULL)) ;
336 }
337 }
338 }
339
340 if (nan_test != NULL)
341 {
342 // use the new operator instead of the built-in one
343 op = nan_test ;
344 b2 = NULL ;
345 }
346 }
347
348 //--------------------------------------------------------------------------
349 // compute C<M> += select (A, b2)
350 //--------------------------------------------------------------------------
351
352 OK1 (C, GxB_Matrix_select (C, M, accum, op, A, (GxB_Scalar) b2, desc)) ;
353
354 //--------------------------------------------------------------------------
355 // free shallow copies
356 //--------------------------------------------------------------------------
357
358 OK (GrB_Matrix_free (&M)) ;
359 OK (GrB_Matrix_free (&A)) ;
360 OK (GrB_Matrix_free (&b)) ;
361 OK (GrB_Descriptor_free (&desc)) ;
362 OK (GrB_BinaryOp_free (&nan_test)) ;
363
364 //--------------------------------------------------------------------------
365 // export the output matrix C back to MATLAB
366 //--------------------------------------------------------------------------
367
368 pargout [0] = gb_export (&C, kind) ;
369 pargout [1] = mxCreateDoubleScalar (kind) ;
370 GB_WRAPUP ;
371 }
372
373