1 //------------------------------------------------------------------------------
2 // GB_kron: C<M> = accum (C, kron(A,B))
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 // C<M> = accum (C, kron(A,B))
11
12 // The input matrices A and B are optionally transposed.
13
14 #include "GB_kron.h"
15 #include "GB_mxm.h"
16 #include "GB_transpose.h"
17 #include "GB_accum_mask.h"
18
19 #define GB_FREE_WORK \
20 { \
21 GB_phbix_free (AT) ; \
22 GB_phbix_free (BT) ; \
23 }
24
25 #define GB_FREE_ALL \
26 { \
27 GB_FREE_WORK ; \
28 GB_phbix_free (T) ; \
29 }
30
GB_kron(GrB_Matrix C,const bool C_replace,const GrB_Matrix M,const bool Mask_comp,const bool Mask_struct,const GrB_BinaryOp accum,const GrB_BinaryOp op_in,const GrB_Matrix A,bool A_transpose,const GrB_Matrix B,bool B_transpose,GB_Context Context)31 GrB_Info GB_kron // C<M> = accum (C, kron(A,B))
32 (
33 GrB_Matrix C, // input/output matrix for results
34 const bool C_replace, // if true, clear C before writing to it
35 const GrB_Matrix M, // optional mask for C, unused if NULL
36 const bool Mask_comp, // if true, use !M
37 const bool Mask_struct, // if true, use the only structure of M
38 const GrB_BinaryOp accum, // optional accum for Z=accum(C,T)
39 const GrB_BinaryOp op_in, // defines '*' for kron(A,B)
40 const GrB_Matrix A, // input matrix
41 bool A_transpose, // if true, use A' instead of A
42 const GrB_Matrix B, // input matrix
43 bool B_transpose, // if true, use B' instead of B
44 GB_Context Context
45 )
46 {
47
48 //--------------------------------------------------------------------------
49 // check inputs
50 //--------------------------------------------------------------------------
51
52 // C may be aliased with M, A, and/or B
53
54 GrB_Info info ;
55 struct GB_Matrix_opaque T_header, AT_header, BT_header ;
56 GrB_Matrix T = GB_clear_static_header (&T_header) ;
57 GrB_Matrix AT = GB_clear_static_header (&AT_header) ;
58 GrB_Matrix BT = GB_clear_static_header (&BT_header) ;
59 GrB_BinaryOp op = op_in ;
60
61 GB_RETURN_IF_NULL_OR_FAULTY (C) ;
62 GB_RETURN_IF_NULL_OR_FAULTY (A) ;
63 GB_RETURN_IF_NULL_OR_FAULTY (B) ;
64 GB_RETURN_IF_FAULTY (M) ;
65 GB_RETURN_IF_NULL_OR_FAULTY (op) ;
66 GB_RETURN_IF_FAULTY_OR_POSITIONAL (accum) ;
67
68 ASSERT_MATRIX_OK (C, "C input for GB_kron", GB0) ;
69 ASSERT_MATRIX_OK_OR_NULL (M, "M for GB_kron", GB0) ;
70 ASSERT_BINARYOP_OK_OR_NULL (accum, "accum for GB_kron", GB0) ;
71 ASSERT_BINARYOP_OK (op, "op for GB_kron", GB0) ;
72 ASSERT_MATRIX_OK (A, "A for GB_kron", GB0) ;
73 ASSERT_MATRIX_OK (B, "B for GB_kron", GB0) ;
74
75 // check domains and dimensions for C<M> = accum (C,T)
76 GB_OK (GB_compatible (C->type, C, M, Mask_struct, accum, op->ztype,
77 Context)) ;
78
79 // T=op(A,B) via op operator, so A and B must be compatible with z=op(a,b)
80 GB_OK (GB_BinaryOp_compatible (op, NULL, A->type, B->type, GB_ignore_code,
81 Context)) ;
82
83 // delete any lingering zombies and assemble any pending tuples in A and B,
84 // so that cnz = nnz(A) * nnz(B) can be computed. Updates of C and M are
85 // done after this check.
86 GB_MATRIX_WAIT (A) ;
87 GB_MATRIX_WAIT (B) ;
88
89 // check the dimensions of C
90 int64_t anrows = (A_transpose) ? GB_NCOLS (A) : GB_NROWS (A) ;
91 int64_t ancols = (A_transpose) ? GB_NROWS (A) : GB_NCOLS (A) ;
92 int64_t bnrows = (B_transpose) ? GB_NCOLS (B) : GB_NROWS (B) ;
93 int64_t bncols = (B_transpose) ? GB_NROWS (B) : GB_NCOLS (B) ;
94 GrB_Index cnrows, cncols, cnz = 0 ;
95 bool ok = GB_Index_multiply (&cnrows, anrows, bnrows) ;
96 ok = ok && GB_Index_multiply (&cncols, ancols, bncols) ;
97 ok = ok && GB_Index_multiply (&cnz, GB_NNZ (A), GB_NNZ (B)) ;
98 if (!ok || GB_NROWS (C) != cnrows || GB_NCOLS (C) != cncols)
99 {
100 GB_ERROR (GrB_DIMENSION_MISMATCH, "%s:\n"
101 "output is " GBd "-by-" GBd "; must be " GBu "-by-" GBu "\n"
102 "first input is " GBd "-by-" GBd "%s with " GBd " entries\n"
103 "second input is " GBd "-by-" GBd "%s with " GBd " entries",
104 ok ? "Dimensions not compatible:" : "Problem too large:",
105 GB_NROWS (C), GB_NCOLS (C), cnrows, cncols,
106 anrows, ancols, A_transpose ? " (transposed)" : "", GB_NNZ (A),
107 bnrows, bncols, B_transpose ? " (transposed)" : "", GB_NNZ (B)) ;
108 }
109
110 // quick return if an empty mask is complemented
111 GB_RETURN_IF_QUICK_MASK (C, C_replace, M, Mask_comp, Mask_struct) ;
112
113 //--------------------------------------------------------------------------
114 // transpose A and B if requested
115 //--------------------------------------------------------------------------
116
117 bool T_is_csc = C->is_csc ;
118 if (T_is_csc != A->is_csc)
119 {
120 // Flip the sense of A_transpose
121 A_transpose = !A_transpose ;
122 }
123 if (T_is_csc != B->is_csc)
124 {
125 // Flip the sense of B_transpose
126 B_transpose = !B_transpose ;
127 }
128
129 if (!T_is_csc)
130 {
131 if (GB_OP_IS_POSITIONAL (op))
132 {
133 // positional ops must be flipped, with i and j swapped
134 op = GB_positional_binop_ijflip (op) ;
135 }
136 }
137
138 // TODO: if A, B are pattern: do not compute values of AT=A', BT=B'
139 bool A_is_pattern, B_is_pattern ;
140 GB_AxB_pattern (&A_is_pattern, &B_is_pattern, false, op->opcode) ;
141
142 if (A_transpose)
143 {
144 // AT = A' and typecast to op->xtype
145 // transpose: typecast, no op, not in-place
146 GBURBLE ("(A transpose) ") ;
147 GB_OK (GB_transpose (&AT, // AT static
148 A_is_pattern ? A->type : op->xtype, T_is_csc,
149 A, NULL, NULL, NULL, false, Context)) ;
150 ASSERT_MATRIX_OK (A , "A after AT kron", GB0) ;
151 ASSERT_MATRIX_OK (AT, "AT kron", GB0) ;
152 }
153
154 if (B_transpose)
155 {
156 // BT = B' and typecast to op->ytype
157 // transpose: typecast, no op, not in-place
158 GBURBLE ("(B transpose) ") ;
159 GB_OK (GB_transpose (&BT, // BT static
160 B_is_pattern ? B->type : op->ytype, T_is_csc,
161 B, NULL, NULL, NULL, false, Context)) ;
162 ASSERT_MATRIX_OK (BT, "BT kron", GB0) ;
163 }
164
165 //--------------------------------------------------------------------------
166 // T = kron(A,B)
167 //--------------------------------------------------------------------------
168
169 GB_OK (GB_kroner (T, T_is_csc, op,
170 A_transpose ? AT : A, A_is_pattern,
171 B_transpose ? BT : B, B_is_pattern, Context)) ;
172
173 GB_FREE_WORK ;
174 ASSERT_MATRIX_OK (T, "T = kron(A,B)", GB0) ;
175
176 //--------------------------------------------------------------------------
177 // C<M> = accum (C,T): accumulate the results into C via the mask
178 //--------------------------------------------------------------------------
179
180 return (GB_accum_mask (C, M, NULL, accum, &T, C_replace, Mask_comp,
181 Mask_struct, Context)) ;
182 }
183
184