1function C = GB_spec_kron (C, Mask, accum, mult, A, B, descriptor)
2%GB_SPEC_KRON a MATLAB mimic of GrB_kronecker
3%
4% Usage:
5% C = GB_spec_kron (C, Mask, accum, mult, A, B, descriptor)
6%
7% Computes C<Mask> = accum(C,T), in GraphBLAS notation, where T = kron(A,B),
8% kron(A',B), kron(A,B') or kron(A',B')
9
10% SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2021, All Rights Reserved.
11% SPDX-License-Identifier: Apache-2.0
12
13%-------------------------------------------------------------------------------
14% get inputs
15%-------------------------------------------------------------------------------
16
17if (nargout > 1 || nargin ~= 7)
18    error ('usage: C = GB_spec_kron (C, Mask, accum, mult, A, B, descriptor)') ;
19end
20
21C = GB_spec_matrix (C) ;
22A = GB_spec_matrix (A) ;
23B = GB_spec_matrix (B) ;
24[mult_op xyclass ztype xtype ytype] = GB_spec_operator (mult, C.class) ;
25[C_replace Mask_comp Atrans Btrans Mask_struct] = ...
26    GB_spec_descriptor (descriptor) ;
27Mask = GB_spec_getmask (Mask, Mask_struct) ;
28
29%-------------------------------------------------------------------------------
30% do the work via a clean MATLAB interpretation of the entire GraphBLAS spec
31%-------------------------------------------------------------------------------
32
33% apply the descriptor to A
34if (Atrans)
35    A.matrix = A.matrix.' ;
36    A.pattern = A.pattern' ;
37end
38
39% apply the descriptor to B
40if (Btrans)
41    B.matrix = B.matrix.' ;
42    B.pattern = B.pattern' ;
43end
44
45% T = A.*B, with typecasting
46[anrows, ancols] = size (A.matrix) ;
47[bnrows, bncols] = size (B.matrix) ;
48cnrows = anrows * bnrows ;
49cncols = ancols * bncols ;
50T.matrix  = GB_spec_zeros ([cnrows cncols], ztype) ;
51T.pattern = false (cnrows, cncols) ;
52
53% first cast the inputs into the x,y types of the operator
54A1 = GB_mex_cast (A.matrix, xtype) ;
55B1 = GB_mex_cast (B.matrix, ytype) ;
56
57for ja = 1:ancols
58    for ia = 1:anrows
59        if (A.pattern (ia,ja))
60            for jb = 1:bncols
61                for ib = 1:bnrows
62                    if (B.pattern (ib,jb))
63                        i = (ia-1) * bnrows + ib ;
64                        j = (ja-1) * bncols + jb ;
65                        if (GB_spec_is_positional (mult))
66                            z = GB_spec_binop_positional (mult_op,ia,ja,ib,jb) ;
67                        else
68                            z = GB_spec_op (mult, A1 (ia,ja), B1 (ib,jb)) ;
69                        end
70                        T.matrix (i,j) = z ;
71                        T.pattern (i,j) = true ;
72                    end
73                end
74            end
75        end
76    end
77end
78
79%{
80    % do the values
81    S = GB_spec_zeros ([bnrows bncols], xtype) ;
82    for j = 1:ancols
83        for i = 1:anrows
84            if A.pattern (i,j)
85                S (:,:) = A1 (i,j) ;
86                ci = (i-1) * bnrows + 1 ;
87                cj = (j-1) * bncols + 1 ;
88                p = B.pattern ;
89                K = GB_spec_op (mult, S(p), B1(p)) ;
90                Tblock = GB_spec_zeros ([bnrows bncols], ztype) ;
91                Tblock (p) = K ;
92                T.matrix  (ci:ci+bnrows-1, cj:cj+bncols-1) = Tblock ;
93                T.pattern (ci:ci+bnrows-1, cj:cj+bncols-1) = B.pattern ;
94            end
95        end
96    end
97%}
98
99assert (isequal (ztype, GB_spec_type (T.matrix))) ;
100T.class = ztype ;
101
102% C<Mask> = accum (C,T): apply the accum, then Mask, and return the result
103C = GB_spec_accum_mask (C, Mask, accum, T, C_replace, Mask_comp, 0) ;
104
105
106