1 //------------------------------------------------------------------------------
2 // gblogassign: logical assignment: C(M) = A
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 // gblogassign computes the MATLAB logical indexing expression C(M) = A. The
11 // matrices C and M must be the same size. M is normally logical but it can be
12 // of any type in this mexFunction. M should not have any explicit zeros. A
13 // is a sparse vector of size nnz(M)-by-1. Scalar expansion is not handled.
14 // Use GrB.subassign (C, M, scalar) for that case.
15
16 // Usage:
17
18 // C = gblogassign (C, M, A)
19
20 // This function is the C equivalent of the following m-function:
21
22 /*
23
24 function C = gblogassign (C, M_input, A)
25 % Computing the MATLAB logical indexing expression C(M) = A in GraphBLAS.
26 % A is a sparse vector of size nnz(M)-by-1 (scalar expansion is not
27 % handled). M is normally a sparse logical matrix, either GraphBLAS or
28 % MATLAB, but it can be of any type. C and M have the same size.
29
30 % make sure all matrices are stored by column
31 save = GrB.format ;
32 GrB.format ('by col') ;
33 M = GrB (m, n, 'logical') ;
34 M = GrB.select (M, '2nd', 'nonzero', M_input) ;
35 if (isequal (GrB.format (A), 'by row'))
36 A = GrB (A) ;
37 end
38
39 [m n] = size (C) ;
40 mnz = nnz (M) ; % A must be mnz-by-1
41 if (~isequal (size (A), [mnz 1]))
42 error ('A must be nnz(M)-by-1')
43 end
44
45 [ai, ~, ax] = GrB.extracttuples (A) ;
46 [mi, mj, ~] = GrB.extracttuples (M) ;
47
48 % construct a subset of the entries of the mask M corresponding to the
49 % entries in A
50 si = mi (ai) ;
51 sj = mj (ai) ;
52 S = GrB.build (si, sj, ax, m, n) ;
53
54 GrB.format (save) ;
55
56 % C<M> = S
57 C = GrB.subassign (C, M, S) ;
58
59 */
60
61
62 // This C mexFunction is faster than the above m-function, since it avoids the
63 // use of GrB.extracttuples. Instead, it accesses the internal structure of the
64 // GrB_Matrix objects. The m-file above is useful for understanding that this
65 // C mexFunction does.
66
67 // C is always returned as a GrB matrix.
68
69 #include "gb_matlab.h"
70
71 #define USAGE "usage: C = gblogassign (C, M, A)"
72 #define ERR "A must be a vector of length nnz(M) for logical indexing, C(M)=A"
73
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])74 void mexFunction
75 (
76 int nargout,
77 mxArray *pargout [ ],
78 int nargin,
79 const mxArray *pargin [ ]
80 )
81 {
82
83 //--------------------------------------------------------------------------
84 // check inputs
85 //--------------------------------------------------------------------------
86
87 gb_usage (nargin == 3 && nargout <= 1, USAGE) ;
88
89 //--------------------------------------------------------------------------
90 // get a deep copy of C, of any sparsity structure
91 //--------------------------------------------------------------------------
92
93 GrB_Matrix C = gb_get_deep (pargin [0]) ;
94 GrB_Index nrows, ncols ;
95 OK (GrB_Matrix_nrows (&nrows, C)) ;
96 OK (GrB_Matrix_ncols (&ncols, C)) ;
97
98 //--------------------------------------------------------------------------
99 // get M
100 //--------------------------------------------------------------------------
101
102 // make M boolean, sparse/hyper, stored by column, and drop explicit zeros
103 GrB_Matrix M_input = gb_get_shallow (pargin [1]) ;
104 GrB_Matrix M = gb_new (GrB_BOOL, nrows, ncols, GxB_BY_COL,
105 GxB_SPARSE + GxB_HYPERSPARSE) ;
106 OK1 (M, GxB_Matrix_select (M, NULL, NULL, GxB_NONZERO, M_input,
107 NULL, NULL)) ;
108
109 OK (GrB_Matrix_free (&M_input)) ;
110 GrB_Index mnz ;
111 OK (GrB_Matrix_nvals (&mnz, M)) ;
112
113 //--------------------------------------------------------------------------
114 // get A
115 //--------------------------------------------------------------------------
116
117 GrB_Matrix A_input = gb_get_shallow (pargin [2]) ;
118 GrB_Matrix A = A_input ;
119 GrB_Type atype ;
120 GrB_Index anrows, ancols, anz ;
121 GxB_Format_Value fmt ;
122 int A_sparsity ;
123 OK (GrB_Matrix_nrows (&anrows, A)) ;
124 OK (GrB_Matrix_ncols (&ancols, A)) ;
125 OK (GxB_Matrix_type (&atype, A)) ;
126 OK (GrB_Matrix_nvals (&anz, A)) ;
127 OK (GxB_Matrix_Option_get (A, GxB_FORMAT, &fmt)) ;
128 OK (GxB_Matrix_Option_get (A, GxB_SPARSITY_STATUS, &A_sparsity)) ;
129
130 GrB_Matrix A_copy = NULL ;
131 GrB_Matrix A_copy2 = NULL ;
132
133 // make sure A is not bitmap; it can be sparse, hypersparse, or full
134 if (A_sparsity == GxB_BITMAP)
135 {
136 OK (GrB_Matrix_dup (&A_copy2, A)) ;
137 OK1 (A_copy2, GxB_Matrix_Option_set (A_copy2, GxB_SPARSITY_CONTROL,
138 GxB_SPARSE + GxB_HYPERSPARSE + GxB_FULL)) ;
139 A = A_copy2 ;
140 }
141
142 // make sure A is a vector of the right size
143 if (mnz == 0)
144 {
145 // M is empty, so A must have no entries. The dimensions and format of
146 // A are not relevant, since the content of A will not be accessed.
147 CHECK_ERROR (anz != 0, ERR) ;
148 }
149 else if (anrows == 1)
150 {
151 // A is 1-by-ancols; ensure it is has length nnz(M), and held by row,
152 // or transpose to ancols-by-1 and held by column.
153 CHECK_ERROR (ancols != mnz, ERR) ;
154 if (fmt == GxB_BY_COL)
155 {
156 // A is 1-by-ancols and held by column: transpose it
157 A_copy = gb_new (atype, mnz, 1, GxB_BY_COL,
158 GxB_SPARSE + GxB_HYPERSPARSE + GxB_FULL) ;
159 OK1 (A_copy, GrB_transpose (A_copy, NULL, NULL, A, NULL)) ;
160 OK1 (A_copy, GrB_Matrix_wait (&A_copy)) ;
161 A = A_copy ;
162 }
163 }
164 else if (ancols == 1)
165 {
166 // A is anrows-by-1; ensure it is has length nnz(M), and held by col
167 // or transpose to 1-by-anrows and held by row.
168 CHECK_ERROR (anrows != mnz, ERR) ;
169 if (fmt == GxB_BY_ROW)
170 {
171 // A is anrows-by-1 and held by row: transpose it
172 A_copy = gb_new (atype, 1, mnz, GxB_BY_ROW,
173 GxB_SPARSE + GxB_HYPERSPARSE + GxB_FULL) ;
174 OK1 (A_copy, GrB_transpose (A_copy, NULL, NULL, A, NULL)) ;
175 OK1 (A_copy, GrB_Matrix_wait (&A_copy)) ;
176 A = A_copy ;
177 }
178 }
179 else
180 {
181 ERROR (ERR) ;
182 }
183
184 //--------------------------------------------------------------------------
185 // extract the values and pattern of A
186 //--------------------------------------------------------------------------
187
188 // TODO: use a shallow variant of GxB*export to access content of M and A
189 GrB_Index *Ai = (GrB_Index *) A->i ;
190 void *Ax = A->x ;
191 double empty ;
192 if (Ax == NULL) Ax = &empty ;
193
194 //--------------------------------------------------------------------------
195 // extract the pattern of M
196 //--------------------------------------------------------------------------
197
198 GrB_Index *Mi = (GrB_Index *) (M->i) ;
199 GrB_Index *Mj = mxMalloc (MAX (mnz, 1) * sizeof (GrB_Index)) ;
200 OK (GrB_Matrix_extractTuples_BOOL (NULL, Mj, NULL, &mnz, M)) ;
201
202 //--------------------------------------------------------------------------
203 // construct a subset of the pattern of M corresponding to the entries of A
204 //--------------------------------------------------------------------------
205
206 GrB_Index *Si = mxMalloc (MAX (anz, 1) * sizeof (GrB_Index)) ;
207 GrB_Index *Sj = mxMalloc (MAX (anz, 1) * sizeof (GrB_Index)) ;
208
209 GB_matlab_helper5 (Si, Sj, Mi, Mj, M->vlen, Ai, A->vlen, anz) ;
210
211 GrB_Matrix S = gb_new (atype, nrows, ncols, GxB_BY_COL, 0) ;
212
213 if (atype == GrB_BOOL)
214 {
215 OK1 (S, GrB_Matrix_build_BOOL (S, Si, Sj, Ax, anz, GrB_LOR)) ;
216 }
217 else if (atype == GrB_INT8)
218 {
219 OK1 (S, GrB_Matrix_build_INT8 (S, Si, Sj, Ax, anz, GrB_PLUS_INT8)) ;
220 }
221 else if (atype == GrB_INT16)
222 {
223 OK1 (S, GrB_Matrix_build_INT16 (S, Si, Sj, Ax, anz, GrB_PLUS_INT16)) ;
224 }
225 else if (atype == GrB_INT32)
226 {
227 OK1 (S, GrB_Matrix_build_INT32 (S, Si, Sj, Ax, anz, GrB_PLUS_INT32)) ;
228 }
229 else if (atype == GrB_INT64)
230 {
231 OK1 (S, GrB_Matrix_build_INT64 (S, Si, Sj, Ax, anz, GrB_PLUS_INT64)) ;
232 }
233 else if (atype == GrB_UINT8)
234 {
235 OK1 (S, GrB_Matrix_build_UINT8 (S, Si, Sj, Ax, anz, GrB_PLUS_UINT8)) ;
236 }
237 else if (atype == GrB_UINT16)
238 {
239 OK1 (S, GrB_Matrix_build_UINT16 (S, Si, Sj, Ax, anz, GrB_PLUS_UINT16)) ;
240 }
241 else if (atype == GrB_UINT32)
242 {
243 OK1 (S, GrB_Matrix_build_UINT32 (S, Si, Sj, Ax, anz, GrB_PLUS_UINT32)) ;
244 }
245 else if (atype == GrB_UINT64)
246 {
247 OK1 (S, GrB_Matrix_build_UINT64 (S, Si, Sj, Ax, anz, GrB_PLUS_UINT64)) ;
248 }
249 else if (atype == GrB_FP32)
250 {
251 OK1 (S, GrB_Matrix_build_FP32 (S, Si, Sj, Ax, anz, GrB_PLUS_FP32)) ;
252 }
253 else if (atype == GrB_FP64)
254 {
255 OK1 (S, GrB_Matrix_build_FP64 (S, Si, Sj, Ax, anz, GrB_PLUS_FP64)) ;
256 }
257 else if (atype == GxB_FC32)
258 {
259 OK1 (S, GxB_Matrix_build_FC32 (S, Si, Sj, Ax, anz, GxB_PLUS_FC32)) ;
260 }
261 else if (atype == GxB_FC64)
262 {
263 OK1 (S, GxB_Matrix_build_FC64 (S, Si, Sj, Ax, anz, GxB_PLUS_FC64)) ;
264 }
265 else
266 {
267 ERROR ("unsupported type") ;
268 }
269
270 OK (GrB_Matrix_free (&A_copy)) ;
271 OK (GrB_Matrix_free (&A_copy2)) ;
272
273 //--------------------------------------------------------------------------
274 // C<M> = S
275 //--------------------------------------------------------------------------
276
277 OK1 (C, GxB_Matrix_subassign (C, M, NULL,
278 S, GrB_ALL, nrows, GrB_ALL, ncols, NULL)) ;
279
280 //--------------------------------------------------------------------------
281 // free shallow copies and temporary matrices
282 //--------------------------------------------------------------------------
283
284 // OK: Si, Sj, and Mj were allocated above from mxMalloc, never in a
285 // GrB_Matrix
286 gb_mxfree (&Si) ;
287 gb_mxfree (&Sj) ;
288 gb_mxfree (&Mj) ;
289 OK (GrB_Matrix_free (&S)) ;
290 OK (GrB_Matrix_free (&M)) ;
291 OK (GrB_Matrix_free (&A_input)) ;
292
293 //--------------------------------------------------------------------------
294 // export the output matrix C back to MATLAB as a GraphBLAS matrix
295 //--------------------------------------------------------------------------
296
297 pargout [0] = gb_export (&C, KIND_GRB) ;
298 GB_WRAPUP ;
299 }
300
301