1 //------------------------------------------------------------------------------
2 // GB_mex_subassign: C(I,J)<M> = accum (C (I,J), A)
3 //------------------------------------------------------------------------------
4
5 // SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2021, All Rights Reserved.
6 // SPDX-License-Identifier: Apache-2.0
7
8 // This function is a wrapper for all GxB_*_subassign functions.
9 // For these uses, the mask M must always be the same size as C(I,J) and A.
10
11 // GxB_Matrix_subassign: M has the same size as C(I,J) and A
12 // GxB_Matrix_subassign_TYPE: M has the same size as C(I,J). A is scalar.
13
14 // GxB_Vector_subassign: M has the same size as C(I,J) and A
15 // GxB_Vector_subassign_TYPE: M has the same size as C(I,J). A is scalar.
16
17 // GxB_Col_subassign: on input to GB_mex_subassign, M and A are a single
18 // columns, the same size as the single subcolumn C(I,j). They are not column
19 // vectors.
20
21 // GxB_Row_subassign: on input to GB_mex_subassign, M and A are single ROWS,
22 // the same size as the single subrow C(i,J). They are not column vectors.
23 // Before GxB_Row_subassign is called, the A and the mask M (if present) are
24 // transposed.
25
26 // Thus, in all cases, A and the mask M (if present) have the same size as
27 // C(I,J), except in the case where A is a scalar. In that case, A is
28 // implicitly expanded into a matrix the same size as C(I,J), but this occurs
29 // inside GxB_*subassign, not here.
30
31 // This function does the same thing as the MATLAB mimic GB_spec_subassign.m.
32
33 //------------------------------------------------------------------------------
34
35 #include "GB_mex.h"
36
37 #define USAGE "[C,s,t] = GB_mex_subassign " \
38 "(C, M, accum, A, I, J, desc, reduce) or (C, Work)"
39
40 #define FREE_ALL \
41 { \
42 bool A_is_M = (A == M) ; \
43 bool A_is_C = (A == C) ; \
44 bool C_is_M = (C == M) ; \
45 GrB_Matrix_free_(&A) ; \
46 if (A_is_C) C = NULL ; \
47 if (A_is_M) M = NULL ; \
48 GrB_Matrix_free_(&C) ; \
49 if (C_is_M) M = NULL ; \
50 GrB_Matrix_free_(&M) ; \
51 GrB_Descriptor_free_(&desc) ; \
52 if (!user_complex) GrB_Monoid_free_(&reduce) ; \
53 GB_mx_put_global (true) ; \
54 }
55
56 #define GET_DEEP_COPY \
57 { \
58 C = GB_mx_mxArray_to_Matrix (pargin [0], "C input", true, true) ; \
59 if (have_sparsity_control) \
60 { \
61 GxB_Matrix_Option_set (C, GxB_SPARSITY_CONTROL, C_sparsity_control) ; \
62 } \
63 if (nargin > 3 && mxIsChar (pargin [1])) \
64 { \
65 M = GB_mx_alias ("M", pargin [1], "C", C, "A", A) ; \
66 } \
67 if (nargin > 3 && mxIsChar (pargin [3])) \
68 { \
69 A = GB_mx_alias ("A", pargin [3], "C", C, "M", M) ; \
70 } \
71 }
72
73 #define FREE_DEEP_COPY \
74 { \
75 if (A == C) A = NULL ; \
76 if (M == C) M = NULL ; \
77 GrB_Matrix_free_(&C) ; \
78 }
79
80 GrB_Matrix C = NULL ;
81 GrB_Matrix M = NULL ;
82 GrB_Matrix A = NULL ;
83 GrB_Matrix mask = NULL, u = NULL ;
84 GrB_Descriptor desc = NULL ;
85 GrB_BinaryOp accum = NULL ;
86 GrB_Index *I = NULL, ni = 0, I_range [3] ;
87 GrB_Index *J = NULL, nj = 0, J_range [3] ;
88 bool ignore ;
89 bool malloc_debug = false ;
90 GrB_Info info = GrB_SUCCESS ;
91 GrB_Monoid reduce = NULL ;
92 GrB_BinaryOp op = NULL ;
93 bool user_complex = false ;
94 int C_sparsity_control ;
95 int M_sparsity_control ;
96 bool have_sparsity_control = false ;
97
98 GrB_Info assign (GB_Context Context) ;
99
100 GrB_Info many_subassign
101 (
102 int nwork,
103 int fA,
104 int fI,
105 int fJ,
106 int faccum,
107 int fM,
108 int fdesc,
109 GrB_Type ctype,
110 const mxArray *pargin [ ],
111 GB_Context Context
112 ) ;
113
114 //------------------------------------------------------------------------------
115 // assign: perform a single assignment
116 //------------------------------------------------------------------------------
117
118 #define OK(method) \
119 { \
120 info = method ; \
121 if (info != GrB_SUCCESS) \
122 { \
123 GrB_Matrix_free_(&mask) ; \
124 GrB_Matrix_free_(&u) ; \
125 return (info) ; \
126 } \
127 }
128
assign(GB_Context Context)129 GrB_Info assign (GB_Context Context)
130 {
131 bool at = (desc != NULL && desc->in0 == GrB_TRAN) ;
132 GrB_Info info ;
133
134 int pr = GB0 ;
135 bool ph = (pr > 0) ;
136
137 ASSERT_MATRIX_OK (C, "C before mex assign", pr) ;
138 ASSERT_BINARYOP_OK_OR_NULL (accum, "accum for mex assign", pr) ;
139 ASSERT_MATRIX_OK (A, "A for mex assign", pr) ;
140
141 if (GB_NROWS (A) == 1 && GB_NCOLS (A) == 1 && GB_NNZ (A) == 1)
142 {
143 // scalar expansion to matrix or vector
144 GB_void *Ax = A->x ;
145
146 if (ni == 1 && nj == 1 && M == NULL && I != GrB_ALL && J != GrB_ALL
147 && GB_op_is_second (accum, C->type) && A->type->code <= GB_FC64_code
148 && desc == NULL)
149 {
150 if (ph) printf ("setElement\n") ;
151 // test GrB_Matrix_setElement
152 #define ASSIGN(prefix,suffix,type) \
153 { \
154 type x = ((type *) Ax) [0] ; \
155 OK (prefix ## Matrix_setElement ## suffix \
156 (C, x, I [0], J [0])) ; \
157 } break ;
158
159 switch (A->type->code)
160 {
161
162 case GB_BOOL_code : ASSIGN (GrB_, _BOOL, bool) ;
163 case GB_INT8_code : ASSIGN (GrB_, _INT8, int8_t) ;
164 case GB_UINT8_code : ASSIGN (GrB_, _UINT8, uint8_t) ;
165 case GB_INT16_code : ASSIGN (GrB_, _INT16, int16_t) ;
166 case GB_UINT16_code : ASSIGN (GrB_, _UINT16, uint16_t) ;
167 case GB_INT32_code : ASSIGN (GrB_, _INT32, int32_t) ;
168 case GB_UINT32_code : ASSIGN (GrB_, _UINT32, uint32_t) ;
169 case GB_INT64_code : ASSIGN (GrB_, _INT64, int64_t) ;
170 case GB_UINT64_code : ASSIGN (GrB_, _UINT64, uint64_t) ;
171 case GB_FP32_code : ASSIGN (GrB_, _FP32, float) ;
172 case GB_FP64_code : ASSIGN (GrB_, _FP64, double) ;
173 case GB_FC32_code : ASSIGN (GxB_, _FC32, GxB_FC32_t) ;
174 case GB_FC64_code : ASSIGN (GxB_, _FC64, GxB_FC64_t) ;
175 case GB_UDT_code :
176 default:
177 FREE_ALL ;
178 mexErrMsgTxt ("unsupported type") ;
179 }
180 #undef ASSIGN
181
182 }
183 if (GB_VECTOR_OK (C) && GB_VECTOR_OK (M))
184 {
185
186 // test GxB_Vector_subassign_scalar functions
187 if (ph) printf ("scalar assign to vector\n") ;
188 #define ASSIGN(suffix,type) \
189 { \
190 type x = ((type *) Ax) [0] ; \
191 OK (GxB_Vector_subassign ## suffix \
192 ((GrB_Vector) C, (GrB_Vector) M, \
193 accum, x, I, ni, desc)) ; \
194 } break ;
195
196 switch (A->type->code)
197 {
198
199 case GB_BOOL_code : ASSIGN (_BOOL, bool) ;
200 case GB_INT8_code : ASSIGN (_INT8, int8_t) ;
201 case GB_UINT8_code : ASSIGN (_UINT8, uint8_t) ;
202 case GB_INT16_code : ASSIGN (_INT16, int16_t) ;
203 case GB_UINT16_code : ASSIGN (_UINT16, uint16_t) ;
204 case GB_INT32_code : ASSIGN (_INT32, int32_t) ;
205 case GB_UINT32_code : ASSIGN (_UINT32, uint32_t) ;
206 case GB_INT64_code : ASSIGN (_INT64, int64_t) ;
207 case GB_UINT64_code : ASSIGN (_UINT64, uint64_t) ;
208 case GB_FP32_code : ASSIGN (_FP32, float) ;
209 case GB_FP64_code : ASSIGN (_FP64, double) ;
210 case GB_FC32_code : ASSIGN (_FC32, GxB_FC32_t) ;
211 case GB_FC64_code : ASSIGN (_FC64, GxB_FC64_t) ;
212 case GB_UDT_code :
213 {
214 // user-defined Complex type
215 OK (GxB_Vector_subassign_UDT
216 ((GrB_Vector) C, (GrB_Vector) M,
217 accum, Ax, I, ni, desc)) ;
218 }
219 break ;
220 default:
221 FREE_ALL ;
222 mexErrMsgTxt ("unsupported type") ;
223 }
224 #undef ASSIGN
225
226 }
227 else
228 {
229
230 // test Matrix_subassign_scalar functions
231 if (ph) printf ("scalar assign to matrix\n") ;
232 #define ASSIGN(suffix,type) \
233 { \
234 type x = ((type *) Ax) [0] ; \
235 OK (GxB_Matrix_subassign ## suffix \
236 (C, M, accum, x, I, ni, J, nj,desc)) ; \
237 } break ;
238
239 switch (A->type->code)
240 {
241
242 case GB_BOOL_code : ASSIGN (_BOOL, bool) ;
243 case GB_INT8_code : ASSIGN (_INT8, int8_t) ;
244 case GB_UINT8_code : ASSIGN (_UINT8, uint8_t) ;
245 case GB_INT16_code : ASSIGN (_INT16, int16_t) ;
246 case GB_UINT16_code : ASSIGN (_UINT16, uint16_t) ;
247 case GB_INT32_code : ASSIGN (_INT32, int32_t) ;
248 case GB_UINT32_code : ASSIGN (_UINT32, uint32_t) ;
249 case GB_INT64_code : ASSIGN (_INT64, int64_t) ;
250 case GB_UINT64_code : ASSIGN (_UINT64, uint64_t) ;
251 case GB_FP32_code : ASSIGN (_FP32, float) ;
252 case GB_FP64_code : ASSIGN (_FP64, double) ;
253 case GB_FC32_code : ASSIGN (_FC32, GxB_FC32_t) ;
254 case GB_FC64_code : ASSIGN (_FC64, GxB_FC64_t) ;
255 case GB_UDT_code :
256 {
257 // user-defined Complex type
258 OK (GxB_Matrix_subassign_UDT
259 (C, M, accum, Ax, I, ni, J, nj, desc)) ;
260 }
261 break ;
262
263 default:
264 FREE_ALL ;
265 mexErrMsgTxt ("unsupported type") ;
266 }
267 #undef ASSIGN
268
269 }
270 }
271 else if (GB_VECTOR_OK (C) && GB_VECTOR_OK (A) &&
272 (M == NULL || GB_VECTOR_OK (M)) && !at)
273 {
274 // test GxB_Vector_subassign
275 if (ph) printf ("vector assign\n") ;
276 OK (GxB_Vector_subassign_((GrB_Vector) C, (GrB_Vector) M, accum,
277 (GrB_Vector) A, I, ni, desc)) ;
278 }
279 else if (GB_VECTOR_OK (A) && nj == 1 &&
280 (M == NULL || GB_VECTOR_OK (M)) && !at)
281 {
282 // test GxB_Col_subassign
283 if (ph) printf ("col assign\n") ;
284 OK (GxB_Col_subassign_(C, (GrB_Vector) M, accum, (GrB_Vector) A,
285 I, ni, J [0], desc)) ;
286 }
287 else if (A->vlen == 1 && ni == 1 && nj > 0 &&
288 (M == NULL || M->vlen == 1) && !at)
289 {
290 // test GxB_Row_subassign; this is not meant to be efficient,
291 // just for testing
292 if (ph) printf ("row assign\n") ;
293 if (M != NULL)
294 {
295 // mask = M'
296 int64_t mnrows, mncols ;
297 OK (GrB_Matrix_nrows (&mnrows, M)) ;
298 OK (GrB_Matrix_ncols (&mncols, M)) ;
299 OK (GrB_Matrix_new (&mask, M->type, mncols, mnrows)) ;
300 OK (GrB_transpose (mask, NULL, NULL, M, NULL)) ;
301 mask->is_csc = true ;
302 ASSERT (GB_VECTOR_OK (mask)) ;
303 }
304 // u = A'
305 int64_t ancols, anrows ;
306 OK (GrB_Matrix_nrows (&anrows, A)) ;
307 OK (GrB_Matrix_ncols (&ancols, A)) ;
308 OK (GrB_Matrix_new (&u, A->type, ancols, anrows)) ;
309 OK (GrB_transpose (u, NULL, NULL, A, NULL)) ;
310 u->is_csc = true ;
311 ASSERT (GB_VECTOR_OK (u)) ;
312 OK (GxB_Row_subassign_(C, (GrB_Vector) mask, accum, (GrB_Vector) u,
313 I [0], J, nj, desc)) ;
314 GrB_Matrix_free_(&mask) ;
315 GrB_Matrix_free_(&u) ;
316 }
317 else
318 {
319 // standard submatrix assignment
320 if (ph) printf ("submatrix assign\n") ;
321 OK (GxB_Matrix_subassign_(C, M, accum, A, I, ni, J, nj, desc)) ;
322 }
323
324 ASSERT_MATRIX_OK (C, "C after assign", pr) ;
325 return (info) ;
326 }
327
328 //------------------------------------------------------------------------------
329 // many_subassign: do a sequence of assignments
330 //------------------------------------------------------------------------------
331
332 // The list of assignments is in a struct array
333
many_subassign(int nwork,int fA,int fI,int fJ,int faccum,int fM,int fdesc,GrB_Type ctype,const mxArray * pargin[],GB_Context Context)334 GrB_Info many_subassign
335 (
336 int nwork,
337 int fA,
338 int fI,
339 int fJ,
340 int faccum,
341 int fM,
342 int fdesc,
343 GrB_Type ctype,
344 const mxArray *pargin [ ],
345 GB_Context Context
346 )
347 {
348 GrB_Info info = GrB_SUCCESS ;
349
350 for (int64_t k = 0 ; k < nwork ; k++)
351 {
352
353 //----------------------------------------------------------------------
354 // get the kth work to do
355 //----------------------------------------------------------------------
356
357 // each struct has fields A, I, J, and optionally Mask, accum, and desc
358
359 mxArray *p ;
360
361 // [ turn off malloc debugging
362 bool save = GB_Global_malloc_debug_get ( ) ;
363 GB_Global_malloc_debug_set (false) ;
364
365 // get M (deep copy)
366 M = NULL ;
367 if (fM >= 0)
368 {
369 p = mxGetFieldByNumber (pargin [1], k, fM) ;
370 M = GB_mx_mxArray_to_Matrix (p, "Mask", true, false) ;
371 if (M == NULL && !mxIsEmpty (p))
372 {
373 FREE_ALL ;
374 mexErrMsgTxt ("M failed") ;
375 }
376 if (have_sparsity_control)
377 {
378 GxB_Matrix_Option_set (M, GxB_SPARSITY_CONTROL,
379 M_sparsity_control) ;
380 }
381 }
382
383 // get A (true copy)
384 p = mxGetFieldByNumber (pargin [1], k, fA) ;
385 A = GB_mx_mxArray_to_Matrix (p, "A", true, true) ;
386 if (A == NULL)
387 {
388 FREE_ALL ;
389 mexErrMsgTxt ("A failed") ;
390 }
391
392 // get accum, if present
393 accum = NULL ;
394 if (faccum >= 0)
395 {
396 p = mxGetFieldByNumber (pargin [1], k, faccum) ;
397 user_complex = (Complex != GxB_FC64)
398 && (C->type == Complex || A->type == Complex) ;
399 if (!GB_mx_mxArray_to_BinaryOp (&accum, p, "accum",
400 C->type, user_complex))
401 {
402 FREE_ALL ;
403 mexErrMsgTxt ("accum failed") ;
404 }
405 }
406
407 // get I
408 p = mxGetFieldByNumber (pargin [1], k, fI) ;
409 if (!GB_mx_mxArray_to_indices (&I, p, &ni, I_range, &ignore))
410 {
411 FREE_ALL ;
412 mexErrMsgTxt ("I failed") ;
413 }
414
415 // get J
416 p = mxGetFieldByNumber (pargin [1], k, fJ) ;
417 if (!GB_mx_mxArray_to_indices (&J, p, &nj, J_range, &ignore))
418 {
419 FREE_ALL ;
420 mexErrMsgTxt ("J failed") ;
421 }
422
423 // get desc
424 desc = NULL ;
425 if (fdesc > 0)
426 {
427 p = mxGetFieldByNumber (pargin [1], k, fdesc) ;
428 if (!GB_mx_mxArray_to_Descriptor (&desc, p, "desc"))
429 {
430 FREE_ALL ;
431 mexErrMsgTxt ("desc failed") ;
432 }
433 }
434
435 // restore malloc debugging to test the method
436 GB_Global_malloc_debug_set (save) ; // ]
437
438 //----------------------------------------------------------------------
439 // C(I,J)<M> = A
440 //----------------------------------------------------------------------
441
442 info = assign (Context) ;
443
444 GrB_Matrix_free_(&A) ;
445 GrB_Matrix_free_(&M) ;
446 GrB_Descriptor_free_(&desc) ;
447
448 if (info != GrB_SUCCESS)
449 {
450 return (info) ;
451 }
452 }
453
454 OK (GrB_Matrix_wait_(&C)) ;
455 return (info) ;
456 }
457
458 //------------------------------------------------------------------------------
459 // GB_mex_subassign mexFunction
460 //------------------------------------------------------------------------------
461
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])462 void mexFunction
463 (
464 int nargout,
465 mxArray *pargout [ ],
466 int nargin,
467 const mxArray *pargin [ ]
468 )
469 {
470
471 C = NULL ;
472 M = NULL ;
473 A = NULL ;
474 mask = NULL ;
475 u = NULL ;
476 desc = NULL ;
477 accum = NULL ;
478 I = NULL ; ni = 0 ;
479 J = NULL ; nj = 0 ;
480 malloc_debug = false ;
481 info = GrB_SUCCESS ;
482 reduce = NULL ;
483 op = NULL ;
484 user_complex = false ;
485 C_sparsity_control = GxB_AUTO_SPARSITY ;
486 M_sparsity_control = GxB_AUTO_SPARSITY ;
487 have_sparsity_control = false ;
488
489 //--------------------------------------------------------------------------
490 // check inputs
491 //--------------------------------------------------------------------------
492
493 malloc_debug = GB_mx_get_global (true) ;
494 A = NULL ;
495 C = NULL ;
496 M = NULL ;
497 desc = NULL ;
498 user_complex = false ;
499 op = NULL ;
500 reduce = NULL ;
501
502 GB_CONTEXT (USAGE) ;
503 if (!((nargout == 1 && (nargin == 2 || nargin == 3 ||
504 nargin == 6 || nargin == 7)) ||
505 ((nargout == 2 || nargout == 3) && nargin == 8)))
506 {
507 mexErrMsgTxt ("Usage: " USAGE) ;
508 }
509
510 if (nargin == 2 || nargin == 3)
511 {
512
513 // get sparsity control if present
514 if (nargin == 3)
515 {
516 int n = mxGetNumberOfElements (pargin [2]) ;
517 if (n != 2) mexErrMsgTxt ("invalid sparsity control") ;
518 have_sparsity_control = true ;
519 double *p = mxGetDoubles (pargin [2]) ;
520 C_sparsity_control = (int) p [0] ;
521 M_sparsity_control = (int) p [1] ;
522 }
523
524 // get C (deep copy)
525 GET_DEEP_COPY ;
526 if (C == NULL)
527 {
528 FREE_ALL ;
529 mexErrMsgTxt ("C failed") ;
530 }
531
532 //----------------------------------------------------------------------
533 // get a list of work to do: a struct array of length nwork
534 //----------------------------------------------------------------------
535
536 // each entry is a struct with fields:
537 // Mask, accum, A, I, J, desc
538
539 if (!mxIsStruct (pargin [1]))
540 {
541 FREE_ALL ;
542 mexErrMsgTxt ("2nd argument must be a struct") ;
543 }
544
545 int nwork = mxGetNumberOfElements (pargin [1]) ;
546 int nf = mxGetNumberOfFields (pargin [1]) ;
547 for (int f = 0 ; f < nf ; f++)
548 {
549 mxArray *p ;
550 for (int k = 0 ; k < nwork ; k++)
551 {
552 p = mxGetFieldByNumber (pargin [1], k, f) ;
553 }
554 }
555
556 int fA = mxGetFieldNumber (pargin [1], "A") ;
557 int fI = mxGetFieldNumber (pargin [1], "I") ;
558 int fJ = mxGetFieldNumber (pargin [1], "J") ;
559 int faccum = mxGetFieldNumber (pargin [1], "accum") ;
560 int fM = mxGetFieldNumber (pargin [1], "Mask") ;
561 int fdesc = mxGetFieldNumber (pargin [1], "desc") ;
562
563 if (fA < 0 || fI < 0 || fJ < 0) mexErrMsgTxt ("A,I,J required") ;
564
565 METHOD (many_subassign (nwork, fA, fI, fJ, faccum, fM, fdesc, C->type,
566 pargin, Context)) ;
567
568 }
569 else
570 {
571
572 //----------------------------------------------------------------------
573 // C(I,J)<M> = A, with a single assignment
574 //----------------------------------------------------------------------
575
576 // get M (deep copy)
577 if (!mxIsChar (pargin [1]))
578 {
579 M = GB_mx_mxArray_to_Matrix (pargin [1], "M", true, false) ;
580 if (M == NULL && !mxIsEmpty (pargin [1]))
581 {
582 FREE_ALL ;
583 mexErrMsgTxt ("M failed") ;
584 }
585 }
586
587 // get A (deep copy)
588 if (!mxIsChar (pargin [3]))
589 {
590 A = GB_mx_mxArray_to_Matrix (pargin [3], "A", true, true) ;
591 if (A == NULL)
592 {
593 FREE_ALL ;
594 mexErrMsgTxt ("A failed") ;
595 }
596 }
597
598 // get C (deep copy)
599 GET_DEEP_COPY ;
600 if (C == NULL)
601 {
602 FREE_ALL ;
603 mexErrMsgTxt ("C failed") ;
604 }
605
606 // get accum, if present
607 user_complex = (Complex != GxB_FC64)
608 && (C->type == Complex || A->type == Complex) ;
609 accum = NULL ;
610 if (!GB_mx_mxArray_to_BinaryOp (&accum, pargin [2], "accum",
611 C->type, user_complex))
612 {
613 FREE_ALL ;
614 mexErrMsgTxt ("accum failed") ;
615 }
616
617 // get I
618 if (!GB_mx_mxArray_to_indices (&I, pargin [4], &ni, I_range, &ignore))
619 {
620 FREE_ALL ;
621 mexErrMsgTxt ("I failed") ;
622 }
623
624 // get J
625 if (!GB_mx_mxArray_to_indices (&J, pargin [5], &nj, J_range, &ignore))
626 {
627 FREE_ALL ;
628 mexErrMsgTxt ("J failed") ;
629 }
630
631 // get desc
632 if (!GB_mx_mxArray_to_Descriptor (&desc, PARGIN (6), "desc"))
633 {
634 FREE_ALL ;
635 mexErrMsgTxt ("desc failed") ;
636 }
637
638 if (nargin == 8 && (nargout == 2 || nargout == 3))
639 {
640 // get reduce operator
641 user_complex = (Complex != GxB_FC64) && (C->type == Complex) ;
642 if (!GB_mx_mxArray_to_BinaryOp (&op, PARGIN (7), "op",
643 C->type, user_complex) || op == NULL)
644 {
645 FREE_ALL ;
646 mexErrMsgTxt ("op failed") ;
647 }
648
649 // get the reduce monoid
650 if (user_complex)
651 {
652 if (op == Complex_plus)
653 {
654 reduce = Complex_plus_monoid ;
655 }
656 else if (op == Complex_times)
657 {
658 reduce = Complex_times_monoid ;
659 }
660 else
661 {
662 FREE_ALL ;
663 mexErrMsgTxt ("user reduce failed") ;
664 }
665 }
666 else
667 {
668 // create the reduce monoid
669 if (!GB_mx_Monoid (&reduce, op, malloc_debug))
670 {
671 FREE_ALL ;
672 mexErrMsgTxt ("reduce failed") ;
673 }
674 }
675 }
676
677 // C(I,J)<M> = A
678 METHOD (assign (Context)) ;
679
680 // apply the reduce monoid
681 if (nargin == 8 && (nargout == 2 || nargout == 3))
682 {
683
684 pargout [1] = GB_mx_create_full (1, 1, C->type) ;
685 GB_void *p = mxGetData (pargout [1]) ;
686
687 #define REDUCE(prefix,suffix,type,zero) \
688 { \
689 type c = zero ; \
690 prefix ## Matrix_reduce ## suffix (&c, NULL, reduce, C, NULL) ;\
691 memcpy (p, &c, sizeof (type)) ; \
692 } \
693 break ;
694
695 double d = 0 ;
696
697 switch (C->type->code)
698 {
699
700 case GB_BOOL_code : REDUCE (GrB_, _BOOL, bool , false);
701 case GB_INT8_code : REDUCE (GrB_, _INT8, int8_t , 0) ;
702 case GB_INT16_code : REDUCE (GrB_, _INT16, int16_t , 0) ;
703 case GB_INT32_code : REDUCE (GrB_, _INT32, int32_t , 0) ;
704 case GB_INT64_code : REDUCE (GrB_, _INT64, int64_t , 0) ;
705 case GB_UINT8_code : REDUCE (GrB_, _UINT8, uint8_t , 0) ;
706 case GB_UINT16_code : REDUCE (GrB_, _UINT16, uint16_t , 0) ;
707 case GB_UINT32_code : REDUCE (GrB_, _UINT32, uint32_t , 0) ;
708 case GB_UINT64_code : REDUCE (GrB_, _UINT64, uint64_t , 0) ;
709 case GB_FP32_code : REDUCE (GrB_, _FP32, float , 0) ;
710 case GB_FP64_code : REDUCE (GrB_, _FP64, double , 0) ;
711 case GB_FC32_code :
712 REDUCE (GxB_, _FC32, GxB_FC32_t, GxB_CMPLXF (0,0)) ;
713 case GB_FC64_code :
714 REDUCE (GxB_, _FC64, GxB_FC64_t, GxB_CMPLX (0,0)) ;
715 case GB_UDT_code :
716 {
717 // user-defined Complex type
718 GxB_FC64_t c = GxB_CMPLX (0,0) ;
719 GrB_Matrix_reduce_UDT_((void *) &c, NULL, reduce,
720 C, NULL) ;
721 memcpy (p, &c, sizeof (GxB_FC64_t)) ;
722 }
723 break ;
724
725 default :
726 FREE_ALL ;
727 mexErrMsgTxt ("unknown type: subassign reduce") ;
728 }
729
730 GrB_Matrix_reduce_FP64_(&d, NULL, GxB_PLUS_FP64_MONOID, C, NULL) ;
731 if (nargout > 2) pargout [2] = mxCreateDoubleScalar (d) ;
732
733 }
734 }
735
736 //--------------------------------------------------------------------------
737 // return C to MATLAB as a struct
738 //--------------------------------------------------------------------------
739
740 ASSERT_MATRIX_OK (C, "Final C before wait", GB0) ;
741 GrB_Matrix_wait_(&C) ;
742
743 if (C == A) A = NULL ; // do not free A if it is aliased to C
744 if (C == M) M = NULL ; // do not free M if it is aliased to C
745 pargout [0] = GB_mx_Matrix_to_mxArray (&C, "C assign result", true) ;
746 FREE_ALL ;
747 }
748
749