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
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
|
//------------------------------------------------------------------------------
// GB_AxB_saxpy_sparsity: determine the sparsity structure for C<M or !M>=A*B
//------------------------------------------------------------------------------
// SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2022, All Rights Reserved.
// SPDX-License-Identifier: Apache-2.0
//------------------------------------------------------------------------------
// Determines the sparsity structure for C, for computing C=A*B, C<M>=A*B, or
// C<!M>=A*B, based on the sparsity structures of C (on input), M, A, and B,
// and whether or not M is complemented.
// TODO: When A or B are bitmapped or full, they can be transposed in-place.
// TODO: give the user control over this decision
//------------------------------------------------------------------------------
#include "GB_AxB_saxpy.h"
void GB_AxB_saxpy_sparsity // determine C_sparsity and method to use
(
// output:
int *C_sparsity, // sparsity structure of C
int *saxpy_method, // saxpy method to use
// input:
const GrB_Matrix M, // optional mask for C, unused if NULL
const bool Mask_comp, // if true, use !M
const GrB_Matrix A, // input A matrix
const GrB_Matrix B, // input B matrix
GB_Context Context
)
{
//--------------------------------------------------------------------------
// determine the sparsity of C
//--------------------------------------------------------------------------
if (B->nvec_nonempty < 0)
{
// B->nvec_nonempty is used to select the method
B->nvec_nonempty = GB_nvec_nonempty (B, Context) ;
}
double bnvec = B->nvec_nonempty ;
double m = (double) A->vlen ;
double n = (double) B->vdim ;
double anz = (double) GB_nnz_held (A) ;
// double bnz = (double) GB_nnz_held (B) ;
int M_sparsity = (M == NULL) ? 0 : GB_sparsity (M) ;
int B_sparsity = GB_sparsity (B) ;
int A_sparsity = GB_sparsity (A) ;
bool M_is_hyper = (M_sparsity == GxB_HYPERSPARSE) ;
bool M_is_sparse = (M_sparsity == GxB_SPARSE) ;
if (M != NULL && !Mask_comp && (M_is_hyper || M_is_sparse))
{
//-----------------------------------------------------
// C <M>= A * B
//-----------------------------------------------------
// hyper sparse any hyper
// hyper hyper any hyper
// sparse hyper any sparse/bitmap/full
// sparse sparse any sparse/bitmap/full
// The non-empty columns of C are a subset of the non-empty columns of
// B, so in general, if B is hypersparse, so is C. If B is sparse,
// bitmap, or full, then C must be sparse, regardless of the sparsity
// of A and B. This is a restriction of GB_AxB_saxpy3.c.
if (B_sparsity == GxB_HYPERSPARSE)
{
(*C_sparsity) = GxB_HYPERSPARSE ;
}
else
{
(*C_sparsity) = GxB_SPARSE ;
}
(*saxpy_method) = GB_SAXPY_METHOD_3 ;
}
else
{
//-----------------------------------------------------
// C = A * B
//-----------------------------------------------------
// hyper . hyper hyper
// hyper . sparse hyper
// hyper . bitmap hyper
// hyper . full hyper
// sparse . hyper sparse
// sparse . sparse sparse
// sparse/bitmap . bitmap sparse
// sparse/bitmap . full sparse
// sparse/bitmap . hyper bitmap
// sparse/bitmap . sparse bitmap
// bitmap . bitmap bitmap
// bitmap . full bitmap
// sparse/bitmap . hyper full
// sparse/bitmap . sparse full
// bitmap . bitmap full
// bitmap (***) . full full
// (***): future, compute C as full
//-----------------------------------------------------
// C <M>= A * B
//-----------------------------------------------------
// hyper any hyper hyper
// hyper any sparse hyper
// hyper any bitmap hyper
// hyper any full hyper
// sparse any hyper sparse
// sparse any sparse sparse
// sparse/bitmap any bitmap sparse
// sparse/bitmap any full sparse
// sparse/bitmap any hyper bitmap
// sparse/bitmap any sparse bitmap
// bitmap any bitmap bitmap
// bitmap any full bitmap
// sparse/bitmap bitmap/full hyper full (*)
// sparse/bitmap bitmap/full sparse full (*)
// bitmap bitmap/full bitmap full (*)
// bitmap bitmap/full full full (*)
// (*): if M hyper/sparse, then C is hyper/sparse; see above
//-----------------------------------------------------
// C <!M>= A * B
//-----------------------------------------------------
// hyper any hyper hyper
// hyper any sparse hyper
// hyper any bitmap hyper
// hyper any full hyper
// sparse any hyper sparse
// sparse any sparse sparse
// sparse/bitmap any bitmap sparse
// sparse/bitmap any full sparse
// sparse/bitmap any hyper bitmap
// sparse/bitmap any sparse bitmap
// bitmap any bitmap bitmap
// bitmap any full bitmap
// sparse/bitmap any hyper full
// sparse/bitmap any sparse full
// bitmap any bitmap full
// bitmap any full full
// If M is complemented, or not complemented and bitmap/full, then C
// has the same sparsity as listed above, except when A and B are both
// full.
// For the cases where C is labelled as hyper/bitmap or sparse/bitmap:
// If m*n is much larger than nnz(A)+nnz(B), then always construct C as
// sparse/hyper, not bitmap. TODO: give the user control over this
// decision.
// TODO: for bitmap*hyper and hyper*bitmap, create a hyper_shallow
// version of the hyper matrix (like dot does), and construct C as
// bitmap. Then expand into C into hyper.
switch (B_sparsity)
{
case GxB_HYPERSPARSE :
// H = any * H
(*C_sparsity) = GxB_HYPERSPARSE ;
break ;
case GxB_SPARSE :
switch (A_sparsity)
{
case GxB_HYPERSPARSE :
case GxB_SPARSE :
// S = {S,H} * S : C has the same sparsity as B
(*C_sparsity) = GxB_SPARSE ;
break ;
case GxB_BITMAP :
case GxB_FULL :
// S = {B,F} * S : if B has many empty columns
// B = {B,F} * S : otherwise C is bitmap
(*C_sparsity) = (bnvec < n/4) ? GxB_SPARSE : GxB_BITMAP;
break ;
default: ;
}
break ;
case GxB_BITMAP :
case GxB_FULL :
switch (A_sparsity)
{
case GxB_HYPERSPARSE :
case GxB_SPARSE :
// S = {S,H} * {B,F} : if A is very sparse
// B = {S,H} * {B,F} : otherwise C is bitmap
(*C_sparsity) = (anz < m/20) ? GxB_SPARSE : GxB_BITMAP ;
break ;
case GxB_BITMAP :
case GxB_FULL :
// B = {B,F} * {B,F} : C is bitmap
(*C_sparsity) = GxB_BITMAP ;
break ;
default: ;
}
break ;
default: ;
}
if ((*C_sparsity) == GxB_HYPERSPARSE || (*C_sparsity) == GxB_SPARSE)
{
(*saxpy_method) = GB_SAXPY_METHOD_3 ;
}
else
{
(*saxpy_method) = GB_SAXPY_METHOD_BITMAP ;
}
}
if ((*C_sparsity) == GxB_HYPERSPARSE || (*C_sparsity) == GxB_SPARSE)
{
// If C is sparse or hypersparse, then it will be computed by
// GB_AxB_saxpy3. For this method, if B is hypersparse, C must also be
// hypersparse. Otherwise C must be sparse. This is a requirement of
// GB_AxB_saxpy3, and is also asserted there.
ASSERT ((*C_sparsity) ==
((B_sparsity == GxB_HYPERSPARSE) ? GxB_HYPERSPARSE : GxB_SPARSE)) ;
}
}
|
