1 //------------------------------------------------------------------------------
2 // gbkronecker: sparse matrix Kronecker product
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 // gbkronecker is an interface to GrB_kronecker
11
12 // Usage:
13
14 // C = gbkronecker (op, A, B)
15 // C = gbkronecker (op, A, B, desc)
16 // C = gbkronecker (Cin, accum, op, A, B, desc)
17 // C = gbkronecker (Cin, M, op, A, B, desc)
18 // C = gbkronecker (Cin, M, accum, op, A, B, desc)
19
20 // If Cin is not present then it is implicitly a matrix with no entries, of the
21 // right size (which depends on A, B, and the descriptor).
22
23 #include "gb_matlab.h"
24
25 #define USAGE "usage: C = GrB.kronecker (Cin, M, accum, op, A, B, desc)"
26
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])27 void mexFunction
28 (
29 int nargout,
30 mxArray *pargout [ ],
31 int nargin,
32 const mxArray *pargin [ ]
33 )
34 {
35
36 //--------------------------------------------------------------------------
37 // check inputs
38 //--------------------------------------------------------------------------
39
40 gb_usage (nargin >= 3 && nargin <= 7 && nargout <= 2, USAGE) ;
41
42 //--------------------------------------------------------------------------
43 // find the arguments
44 //--------------------------------------------------------------------------
45
46 mxArray *Matrix [4], *String [2], *Cell [2] ;
47 base_enum_t base ;
48 kind_enum_t kind ;
49 GxB_Format_Value fmt ;
50 int nmatrices, nstrings, ncells, sparsity ;
51 GrB_Descriptor desc ;
52 gb_get_mxargs (nargin, pargin, USAGE, Matrix, &nmatrices, String, &nstrings,
53 Cell, &ncells, &desc, &base, &kind, &fmt, &sparsity) ;
54
55 CHECK_ERROR (nmatrices < 2 || nstrings < 1 || ncells > 0, USAGE) ;
56
57 //--------------------------------------------------------------------------
58 // get the matrices
59 //--------------------------------------------------------------------------
60
61 GrB_Type atype, btype, ctype = NULL ;
62 GrB_Matrix C = NULL, M = NULL, A, B ;
63
64 if (nmatrices == 2)
65 {
66 A = gb_get_shallow (Matrix [0]) ;
67 B = gb_get_shallow (Matrix [1]) ;
68 }
69 else if (nmatrices == 3)
70 {
71 C = gb_get_deep (Matrix [0]) ;
72 A = gb_get_shallow (Matrix [1]) ;
73 B = gb_get_shallow (Matrix [2]) ;
74 }
75 else // if (nmatrices == 4)
76 {
77 C = gb_get_deep (Matrix [0]) ;
78 M = gb_get_shallow (Matrix [1]) ;
79 A = gb_get_shallow (Matrix [2]) ;
80 B = gb_get_shallow (Matrix [3]) ;
81 }
82
83 OK (GxB_Matrix_type (&atype, A)) ;
84 OK (GxB_Matrix_type (&btype, B)) ;
85 if (C != NULL)
86 {
87 OK (GxB_Matrix_type (&ctype, C)) ;
88 }
89
90 //--------------------------------------------------------------------------
91 // get the operators
92 //--------------------------------------------------------------------------
93
94 GrB_BinaryOp accum = NULL, op ;
95
96 if (nstrings == 1)
97 {
98 op = gb_mxstring_to_binop (String [0], atype, btype) ;
99 }
100 else
101 {
102 // if accum appears, then Cin must also appear
103 CHECK_ERROR (C == NULL, USAGE) ;
104 accum = gb_mxstring_to_binop (String [0], ctype, ctype) ;
105 op = gb_mxstring_to_binop (String [1], atype, btype) ;
106 }
107
108 //--------------------------------------------------------------------------
109 // construct C if not present on input
110 //--------------------------------------------------------------------------
111
112 // If C is NULL, then it is not present on input.
113 // Construct C of the right size and type.
114
115 if (C == NULL)
116 {
117 // get the descriptor contents to determine if A and B are transposed
118 GrB_Desc_Value in0, in1 ;
119 OK (GxB_Desc_get (desc, GrB_INP0, &in0)) ;
120 OK (GxB_Desc_get (desc, GrB_INP1, &in1)) ;
121 bool A_transpose = (in0 == GrB_TRAN) ;
122 bool B_transpose = (in1 == GrB_TRAN) ;
123
124 // get the size of A and B
125 GrB_Index anrows, ancols, bnrows, bncols ;
126 if (A_transpose)
127 {
128 OK (GrB_Matrix_nrows (&ancols, A)) ;
129 OK (GrB_Matrix_ncols (&anrows, A)) ;
130 }
131 else
132 {
133 OK (GrB_Matrix_nrows (&anrows, A)) ;
134 OK (GrB_Matrix_ncols (&ancols, A)) ;
135 }
136 if (B_transpose)
137 {
138 OK (GrB_Matrix_nrows (&bncols, B)) ;
139 OK (GrB_Matrix_ncols (&bnrows, B)) ;
140 }
141 else
142 {
143 OK (GrB_Matrix_nrows (&bnrows, B)) ;
144 OK (GrB_Matrix_ncols (&bncols, B)) ;
145 }
146
147 // determine the size of C
148 GrB_Index cnrows = anrows * bnrows ;
149 GrB_Index cncols = ancols * bncols ;
150
151 // use the ztype of the op as the type of C
152 OK (GxB_BinaryOp_ztype (&ctype, op)) ;
153
154 // create the matrix C and set its format and sparsity
155 fmt = gb_get_format (cnrows, cncols, A, B, fmt) ;
156 sparsity = gb_get_sparsity (A, B, sparsity) ;
157 C = gb_new (ctype, cnrows, cncols, fmt, sparsity) ;
158 }
159
160 //--------------------------------------------------------------------------
161 // compute C<M> += kron (A,B)
162 //--------------------------------------------------------------------------
163
164 OK1 (C, GrB_Matrix_kronecker_BinaryOp (C, M, accum, op, A, B, desc)) ;
165
166 //--------------------------------------------------------------------------
167 // free shallow copies
168 //--------------------------------------------------------------------------
169
170 OK (GrB_Matrix_free (&M)) ;
171 OK (GrB_Matrix_free (&A)) ;
172 OK (GrB_Matrix_free (&B)) ;
173 OK (GrB_Descriptor_free (&desc)) ;
174
175 //--------------------------------------------------------------------------
176 // export the output matrix C back to MATLAB
177 //--------------------------------------------------------------------------
178
179 pargout [0] = gb_export (&C, kind) ;
180 pargout [1] = mxCreateDoubleScalar (kind) ;
181 GB_WRAPUP ;
182 }
183
184