1 //------------------------------------------------------------------------------
2 // GB_mx_build_template: build a sparse vector or matrix
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 // This is not a stand-alone function; it is #include'd in
11 // GB_mex_Matrix_build.c and GB_mex_Vector_build.c.
12
13 // The function works the same as C = sparse (I,J,X,nrows,ncols) in MATLAB,
14 // except an optional operator, dup, can be provided. The operator defines
15 // how duplicates are assembled.
16
17 // The type of dup and C may differ. A matrix T is first built that has the
18 // same type as dup, and then typecasted to the desired type of C, given by
19 // the last argument, "type".
20
21 // This particular GraphBLAS implementation provides a well-defined order of
22 // 'summation'. Entries in [I,J,X] are first sorted in increasing order of row
23 // index via a stable sort, with ties broken by the position of the tuple in
24 // the [I,J,X] list. If duplicates appear, they are 'summed' in the order they
25 // appear in the [I,J,X] input. That is, if the same indices i and j appear in
26 // positions k1, k2, k3, and k3 in [I,J,X], where k1 < k2 < k3 < k4, then the
27 // following operations will occur in order:
28
29 // T (i,j) = X (k1) ;
30 // T (i,j) = dup (T (i,j), X (k2)) ;
31 // T (i,j) = dup (T (i,j), X (k3)) ;
32 // T (i,j) = dup (T (i,j), X (k4)) ;
33
34 // Thus, dup need not be associative, and the results are still well-defined.
35 // Using the FIRST operator (a string 'first') means the first value (X (k1))
36 // is used and the rest are ignored. The SECOND operator means the last value
37 // (X (k4)) is used instead.
38
39 #ifdef MATRIX
40 #define MAX_NARGIN 8
41 #define MIN_NARGIN 3
42 #define USAGE "GB_mex_Matrix_build (I,J,X,nrows,ncols,dup,type,csc)"
43 #define I_ARG 0
44 #define J_ARG 1
45 #define X_ARG 2
46 #define NROWS_ARG 3
47 #define NCOLS_ARG 4
48 #define DUP_ARG 5
49 #define TYPE_ARG 6
50 #define CSC_ARG 7
51 #define FREE_WORK GrB_Matrix_free_(&C) ;
52 #else
53 #define MAX_NARGIN 5
54 #define MIN_NARGIN 2
55 #define USAGE "GB_mex_Vector_build (I,X,nrows,dup,type)"
56 #define I_ARG 0
57 #define X_ARG 1
58 #define NROWS_ARG 2
59 #define DUP_ARG 3
60 #define TYPE_ARG 4
61 #define FREE_WORK GrB_Vector_free_(&C) ;
62 #endif
63
64 #define FREE_ALL \
65 { \
66 FREE_WORK ; \
67 GB_mx_put_global (true) ; \
68 }
69
70 #define GET_DEEP_COPY ;
71 #define FREE_DEEP_COPY ;
72
73 bool malloc_debug = false ;
74
75 GrB_Info builder
76 (
77 #ifdef MATRIX
78 GrB_Matrix *Chandle,
79 #else
80 GrB_Vector *Chandle,
81 #endif
82 GrB_Type ctype,
83 GrB_Index nrows,
84 GrB_Index ncols,
85 GrB_Index *I,
86 GrB_Index *J,
87 GB_void *X,
88 GrB_Index ni,
89 GrB_BinaryOp dup,
90 bool C_is_csc,
91 GrB_Type xtype,
92 GB_Context Context
93 ) ;
94
95 //------------------------------------------------------------------------------
96
builder(GrB_Matrix * Chandle,GrB_Type ctype,GrB_Index nrows,GrB_Index ncols,GrB_Index * I,GrB_Index * J,GB_void * X,GrB_Index ni,GrB_BinaryOp dup,bool C_is_csc,GrB_Type xtype,GB_Context Context)97 GrB_Info builder
98 (
99 #ifdef MATRIX
100 GrB_Matrix *Chandle,
101 #else
102 GrB_Vector *Chandle,
103 #endif
104 GrB_Type ctype,
105 GrB_Index nrows,
106 GrB_Index ncols,
107 GrB_Index *I,
108 GrB_Index *J,
109 GB_void *X,
110 GrB_Index ni,
111 GrB_BinaryOp dup,
112 bool C_is_csc,
113 GrB_Type xtype,
114 GB_Context Context
115 )
116 {
117
118 GrB_Info info ;
119 (*Chandle) = NULL ;
120
121 // create the GraphBLAS output object C
122 int sparsity = GxB_SPARSE + GxB_HYPERSPARSE ;
123 #ifdef MATRIX
124 if (C_is_csc)
125 {
126 // create a hypersparse CSC matrix
127 info = GB_new (Chandle, false, // sparse/hyper, new mx header
128 ctype, nrows, ncols, GB_Ap_calloc,
129 true, sparsity, GxB_HYPER_DEFAULT, 1, Context) ;
130 }
131 else
132 {
133 // create a hypersparse CSR matrix
134 info = GB_new (Chandle, false, // sparse/hyper, new mx header
135 ctype, ncols, nrows, GB_Ap_calloc,
136 false, sparsity, GxB_HYPER_DEFAULT, 1, Context) ;
137 }
138 #else
139 info = GrB_Vector_new (Chandle, ctype, nrows) ;
140 #endif
141
142 #ifdef MATRIX
143 GrB_Matrix C = (*Chandle) ;
144 #else
145 GrB_Vector C = (*Chandle) ;
146 #endif
147
148 if (info != GrB_SUCCESS)
149 {
150 FREE_WORK ;
151 return (info) ;
152 }
153
154 // build the matrix or vector from the tuples
155 #ifdef MATRIX
156 #define BUILD(prefix,suffix,type) \
157 info = prefix ## Matrix_build ## suffix (C,I,J,(const type *)X,ni,dup)
158 #else
159 #define BUILD(prefix,suffix,type) \
160 info = prefix ## Vector_build ## suffix (C,I, (const type *)X,ni,dup)
161 #endif
162
163 ASSERT_TYPE_OK (ctype, "ctype for build", GB0) ;
164 ASSERT_BINARYOP_OK (dup, "dup for build", GB0) ;
165 // printf ("code %d biulding ni "GBd"\n", ctype->code, ni) ;
166
167 switch (xtype->code)
168 {
169 case GB_BOOL_code : BUILD (GrB_, _BOOL, bool ) ; break ;
170 case GB_INT8_code : BUILD (GrB_, _INT8, int8_t ) ; break ;
171 case GB_UINT8_code : BUILD (GrB_, _UINT8, uint8_t ) ; break ;
172 case GB_INT16_code : BUILD (GrB_, _INT16, int16_t ) ; break ;
173 case GB_UINT16_code : BUILD (GrB_, _UINT16, uint16_t) ; break ;
174 case GB_INT32_code : BUILD (GrB_, _INT32, int32_t ) ; break ;
175 case GB_UINT32_code : BUILD (GrB_, _UINT32, uint32_t) ; break ;
176 case GB_INT64_code : BUILD (GrB_, _INT64, int64_t ) ; break ;
177 case GB_UINT64_code : BUILD (GrB_, _UINT64, uint64_t) ; break ;
178 case GB_FP32_code : BUILD (GrB_, _FP32, float ) ; break ;
179 case GB_FP64_code : BUILD (GrB_, _FP64, double ) ; break ;
180 case GB_FC32_code : BUILD (GxB_, _FC32, GxB_FC32_t) ; break ;
181 case GB_FC64_code : BUILD (GxB_, _FC64, GxB_FC64_t) ; break ;
182 default :
183 FREE_WORK ;
184 mexErrMsgTxt ("xtype not supported") ;
185 }
186
187 if (info == GrB_SUCCESS)
188 {
189 ASSERT_MATRIX_OK (C, "C built", GB0) ;
190 }
191 else
192 {
193 FREE_WORK ;
194 }
195
196 return (info) ;
197 }
198
199 //------------------------------------------------------------------------------
200
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])201 void mexFunction
202 (
203 int nargout,
204 mxArray *pargout [ ],
205 int nargin,
206 const mxArray *pargin [ ]
207 )
208 {
209
210 GrB_Info info ;
211 malloc_debug = GB_mx_get_global (true) ;
212 GrB_Index *I = NULL, ni = 0, I_range [3] ;
213 GrB_Index *J = NULL, nj = 0, J_range [3] ;
214 bool is_list ;
215 #ifdef MATRIX
216 GrB_Matrix C = NULL ;
217 #else
218 GrB_Vector C = NULL ;
219 #endif
220
221 GB_CONTEXT (USAGE) ;
222
223 // check inputs
224 if (nargout > 1 || nargin < MIN_NARGIN || nargin > MAX_NARGIN)
225 {
226 mexErrMsgTxt ("Usage: C = " USAGE) ;
227 }
228
229 // get I
230 if (!GB_mx_mxArray_to_indices (&I, pargin [I_ARG], &ni, I_range, &is_list))
231 {
232 FREE_ALL ;
233 mexErrMsgTxt ("I failed") ;
234 }
235 if (!is_list)
236 {
237 mexErrMsgTxt ("I is invalid; must be a list") ;
238 }
239
240 #ifdef MATRIX
241 // get J for a matrix
242 if (!GB_mx_mxArray_to_indices (&J, pargin [J_ARG], &nj, J_range, &is_list))
243 {
244 FREE_ALL ;
245 mexErrMsgTxt ("J failed") ;
246 }
247 if (!is_list)
248 {
249 mexErrMsgTxt ("J is invalid; must be a list") ;
250 }
251
252 if (ni != nj)
253 {
254 FREE_ALL ;
255 mexErrMsgTxt ("I and J must be the same size") ;
256 }
257 #endif
258
259 // get X
260 if (ni != mxGetNumberOfElements (pargin [X_ARG]))
261 {
262 FREE_ALL ;
263 mexErrMsgTxt ("I and X must be the same size") ;
264 }
265 if (!(mxIsNumeric (pargin [X_ARG]) || mxIsLogical (pargin [X_ARG])))
266 {
267 FREE_ALL ;
268 mexErrMsgTxt ("X must be a numeric or logical array") ;
269 }
270 if (mxIsSparse (pargin [X_ARG]))
271 {
272 FREE_ALL ;
273 mexErrMsgTxt ("X cannot be sparse") ;
274 }
275 GB_void *X = mxGetData (pargin [X_ARG]) ;
276 GrB_Type xtype = GB_mx_Type (pargin [X_ARG]) ;
277 if (xtype == NULL)
278 {
279 FREE_ALL ;
280 mexErrMsgTxt ("X must be numeric") ;
281 }
282
283 // get the number of rows
284 uint64_t nrows = 0 ;
285 if (nargin > NROWS_ARG)
286 {
287 nrows = (uint64_t) mxGetScalar (pargin [NROWS_ARG]) ;
288 }
289 else
290 {
291 for (int64_t k = 0 ; k < ni ; k++)
292 {
293 nrows = GB_IMAX (nrows, I [k]) ;
294 }
295 nrows++ ;
296 }
297
298 // get the number of columns of a matrix
299 uint64_t ncols = 1 ;
300 #ifdef MATRIX
301 if (nargin > NCOLS_ARG)
302 {
303 ncols = (uint64_t) mxGetScalar (pargin [NCOLS_ARG]) ;
304 }
305 else
306 {
307 ncols = 0 ;
308 for (int64_t k = 0 ; k < ni ; k++)
309 {
310 ncols = GB_IMAX (ncols, J [k]) ;
311 }
312 ncols++ ;
313 }
314 #endif
315
316 // get dup; default is xtype
317 bool user_complex = (Complex != GxB_FC64) && (xtype == Complex) ;
318 GrB_BinaryOp dup ;
319 if (!GB_mx_mxArray_to_BinaryOp (&dup, PARGIN (DUP_ARG), "dup",
320 xtype, user_complex))
321 {
322 FREE_ALL ;
323 mexErrMsgTxt ("dup failed") ;
324 }
325
326 // get the type for C, default is same as xtype
327 GrB_Type ctype = GB_mx_string_to_Type (PARGIN (TYPE_ARG), xtype) ;
328
329 if (dup == NULL)
330 {
331 switch (xtype->code)
332 {
333 case GB_BOOL_code : dup = GrB_PLUS_BOOL ; break ;
334 case GB_INT8_code : dup = GrB_PLUS_INT8 ; break ;
335 case GB_INT16_code : dup = GrB_PLUS_INT16 ; break ;
336 case GB_INT32_code : dup = GrB_PLUS_INT32 ; break ;
337 case GB_INT64_code : dup = GrB_PLUS_INT64 ; break ;
338 case GB_UINT8_code : dup = GrB_PLUS_UINT8 ; break ;
339 case GB_UINT16_code : dup = GrB_PLUS_UINT16 ; break ;
340 case GB_UINT32_code : dup = GrB_PLUS_UINT32 ; break ;
341 case GB_UINT64_code : dup = GrB_PLUS_UINT64 ; break ;
342 case GB_FP32_code : dup = GrB_PLUS_FP32 ; break ;
343 case GB_FP64_code : dup = GrB_PLUS_FP64 ; break ;
344 case GB_FC32_code : dup = GxB_PLUS_FC32 ; break ;
345 case GB_FC64_code : dup = GxB_PLUS_FC64 ; break ;
346 default :
347 mexErrMsgTxt ("unknown operator") ;
348 }
349 }
350
351 bool C_is_csc = true ;
352 #ifdef MATRIX
353 // get the CSC/CSR format
354 if (nargin > CSC_ARG)
355 {
356 C_is_csc = (bool) mxGetScalar (pargin [CSC_ARG]) ;
357 }
358 #endif
359
360 METHOD (builder (&C, ctype, nrows, ncols, I, J, X, ni, dup,
361 C_is_csc, xtype, Context)) ;
362
363 ASSERT_MATRIX_OK (C, "C built", GB0) ;
364
365 // return C to MATLAB as a struct and free the GraphBLAS C
366 #ifdef MATRIX
367 pargout [0] = GB_mx_Matrix_to_mxArray (&C, "C output", true) ;
368 #else
369 pargout [0] = GB_mx_Vector_to_mxArray (&C, "C output", true) ;
370 #endif
371
372 FREE_ALL ;
373 }
374
375