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