1 //------------------------------------------------------------------------------
2 // gbextract: extract entries into 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 // gbextract is an interface to GrB_Matrix_extract and GrB_Matrix_extract_[TYPE],
11 // computing the GraphBLAS expression:
12 
13 //      C<#M,replace> = accum (C, A (I,J)) or
14 //      C<#M,replace> = accum (C, AT (I,J))
15 
16 // Usage:
17 
18 //      C = gbextract (Cin, M, accum, A, I, J, desc)
19 
20 // A is required.  See GrB.m for more details.
21 // If accum or M is used, then Cin must appear.
22 
23 #include "gb_matlab.h"
24 #include "GB_ij.h"
25 
26 #define USAGE "usage: C = GrB.extract (Cin, M, accum, A, I, J, desc)"
27 
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])28 void mexFunction
29 (
30     int nargout,
31     mxArray *pargout [ ],
32     int nargin,
33     const mxArray *pargin [ ]
34 )
35 {
36 
37     //--------------------------------------------------------------------------
38     // check inputs
39     //--------------------------------------------------------------------------
40 
41     gb_usage (nargin >= 1 && nargin <= 7 && nargout <= 2, USAGE) ;
42 
43     //--------------------------------------------------------------------------
44     // find the arguments
45     //--------------------------------------------------------------------------
46 
47     mxArray *Matrix [4], *String [2], *Cell [2] ;
48     base_enum_t base ;
49     kind_enum_t kind ;
50     GxB_Format_Value fmt ;
51     int nmatrices, nstrings, ncells, sparsity ;
52     GrB_Descriptor desc ;
53     gb_get_mxargs (nargin, pargin, USAGE, Matrix, &nmatrices, String, &nstrings,
54         Cell, &ncells, &desc, &base, &kind, &fmt, &sparsity) ;
55 
56     CHECK_ERROR (nmatrices < 1 || nmatrices > 3 || nstrings > 1, USAGE) ;
57 
58     //--------------------------------------------------------------------------
59     // get the matrices
60     //--------------------------------------------------------------------------
61 
62     GrB_Type atype, ctype = NULL ;
63     GrB_Matrix C = NULL, M = NULL, A ;
64 
65     if (nmatrices == 1)
66     {
67         A = gb_get_shallow (Matrix [0]) ;
68     }
69     else if (nmatrices == 2)
70     {
71         C = gb_get_deep    (Matrix [0]) ;
72         A = gb_get_shallow (Matrix [1]) ;
73     }
74     else // if (nmatrices == 3)
75     {
76         C = gb_get_deep    (Matrix [0]) ;
77         M = gb_get_shallow (Matrix [1]) ;
78         A = gb_get_shallow (Matrix [2]) ;
79     }
80 
81     OK (GxB_Matrix_type (&atype, A)) ;
82     if (C != NULL)
83     {
84         OK (GxB_Matrix_type (&ctype, C)) ;
85     }
86 
87     //--------------------------------------------------------------------------
88     // get the operator
89     //--------------------------------------------------------------------------
90 
91     GrB_BinaryOp accum = NULL ;
92 
93     if (nstrings == 1)
94     {
95         // if accum appears, then Cin must also appear
96         CHECK_ERROR (C == NULL, USAGE) ;
97         accum  = gb_mxstring_to_binop  (String [0], ctype, ctype) ;
98     }
99 
100     //--------------------------------------------------------------------------
101     // get the size of A
102     //--------------------------------------------------------------------------
103 
104     GrB_Desc_Value in0 ;
105     OK (GxB_Desc_get (desc, GrB_INP0, &in0)) ;
106     GrB_Index anrows, ancols ;
107     bool A_transpose = (in0 == GrB_TRAN) ;
108     if (A_transpose)
109     {
110         // T = AT (I,J) is to be extracted where AT = A'
111         OK (GrB_Matrix_nrows (&ancols, A)) ;
112         OK (GrB_Matrix_ncols (&anrows, A)) ;
113     }
114     else
115     {
116         // T = A (I,J) is to be extracted
117         OK (GrB_Matrix_nrows (&anrows, A)) ;
118         OK (GrB_Matrix_ncols (&ancols, A)) ;
119     }
120 
121     //--------------------------------------------------------------------------
122     // get I and J
123     //--------------------------------------------------------------------------
124 
125     GrB_Index *I = (GrB_Index *) GrB_ALL ;
126     GrB_Index *J = (GrB_Index *) GrB_ALL ;
127     GrB_Index ni = anrows, nj = ancols ;
128     bool I_allocated = false, J_allocated = false ;
129 
130     if (anrows == 1 && ncells == 1)
131     {
132         // only J is present
133         J = gb_mxcell_to_index (Cell [0], base, ancols, &J_allocated, &nj) ;
134     }
135     else if (ncells == 1)
136     {
137         // only I is present
138         I = gb_mxcell_to_index (Cell [0], base, anrows, &I_allocated, &ni) ;
139     }
140     else if (ncells == 2)
141     {
142         // both I and J are present
143         I = gb_mxcell_to_index (Cell [0], base, anrows, &I_allocated, &ni) ;
144         J = gb_mxcell_to_index (Cell [1], base, ancols, &J_allocated, &nj) ;
145     }
146 
147     //--------------------------------------------------------------------------
148     // construct C if not present on input
149     //--------------------------------------------------------------------------
150 
151     if (C == NULL)
152     {
153         // Cin is not present: determine its size, same type as A.
154         // T = A(I,J) or AT(I,J) will be extracted.
155         // accum must be null
156         int I_kind, J_kind ;
157         int64_t I_colon [3], J_colon [3] ;
158         GrB_Index cnrows, cncols ;
159         GB_ijlength (I, ni, anrows, (int64_t *) &cnrows, &I_kind, I_colon) ;
160         GB_ijlength (J, nj, ancols, (int64_t *) &cncols, &J_kind, J_colon) ;
161         ctype = atype ;
162 
163         // create the matrix C and set its format and sparsity
164         fmt = gb_get_format (cnrows, cncols, A, NULL, fmt) ;
165         sparsity = gb_get_sparsity (A, NULL, sparsity) ;
166         C = gb_new (ctype, cnrows, cncols, fmt, sparsity) ;
167     }
168 
169     //--------------------------------------------------------------------------
170     // compute C<M> += A(I,J) or AT(I,J)
171     //--------------------------------------------------------------------------
172 
173     OK1 (C, GrB_Matrix_extract (C, M, accum, A, I, ni, J, nj, desc)) ;
174 
175     //--------------------------------------------------------------------------
176     // free shallow copies
177     //--------------------------------------------------------------------------
178 
179     OK (GrB_Matrix_free (&M)) ;
180     OK (GrB_Matrix_free (&A)) ;
181     OK (GrB_Descriptor_free (&desc)) ;
182     if (I_allocated) gb_mxfree (&I) ;
183     if (J_allocated) gb_mxfree (&J) ;
184 
185     //--------------------------------------------------------------------------
186     // export the output matrix C back to MATLAB
187     //--------------------------------------------------------------------------
188 
189     pargout [0] = gb_export (&C, kind) ;
190     pargout [1] = mxCreateDoubleScalar (kind) ;
191     GB_WRAPUP ;
192 }
193 
194