1 //------------------------------------------------------------------------------
2 // GB_AxB_rowscale: C = D*B, row scale with diagonal matrix D
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 #include "GB_mxm.h"
11 #include "GB_binop.h"
12 #include "GB_apply.h"
13 #ifndef GBCOMPACT
14 #include "GB_binop__include.h"
15 #endif
16
17 #define GB_FREE_ALL GB_phbix_free (C) ;
18
GB_AxB_rowscale(GrB_Matrix C,const GrB_Matrix D,const GrB_Matrix B,const GrB_Semiring semiring,const bool flipxy,GB_Context Context)19 GrB_Info GB_AxB_rowscale // C = D*B, row scale with diagonal D
20 (
21 GrB_Matrix C, // output matrix, static header
22 const GrB_Matrix D, // diagonal input matrix
23 const GrB_Matrix B, // input matrix
24 const GrB_Semiring semiring, // semiring that defines C=D*A
25 const bool flipxy, // if true, do z=fmult(b,a) vs fmult(a,b)
26 GB_Context Context
27 )
28 {
29
30 //--------------------------------------------------------------------------
31 // check inputs
32 //--------------------------------------------------------------------------
33
34 GrB_Info info ;
35 ASSERT (C != NULL && C->static_header) ;
36
37 ASSERT_MATRIX_OK (D, "D for rowscale A*D", GB0) ;
38 ASSERT (!GB_ZOMBIES (D)) ;
39 ASSERT (!GB_JUMBLED (D)) ;
40 ASSERT (!GB_PENDING (D)) ;
41
42 ASSERT_MATRIX_OK (B, "B for rowscale A*D", GB0) ;
43 ASSERT (!GB_ZOMBIES (B)) ;
44 ASSERT (GB_JUMBLED_OK (B)) ;
45 ASSERT (!GB_PENDING (B)) ;
46
47 ASSERT_SEMIRING_OK (semiring, "semiring for numeric D*A", GB0) ;
48 ASSERT (D->vdim == B->vlen) ;
49 ASSERT (GB_is_diagonal (D, Context)) ;
50
51 ASSERT (!GB_IS_BITMAP (D)) ; // bitmap or full: not needed
52 ASSERT (!GB_IS_BITMAP (B)) ;
53 ASSERT (!GB_IS_FULL (D)) ;
54
55 GBURBLE ("(%s=%s*%s) ",
56 GB_sparsity_char_matrix (B), // C has the sparsity structure of B
57 GB_sparsity_char_matrix (D),
58 GB_sparsity_char_matrix (B)) ;
59
60 //--------------------------------------------------------------------------
61 // get the semiring operators
62 //--------------------------------------------------------------------------
63
64 GrB_BinaryOp mult = semiring->multiply ;
65 ASSERT (mult->ztype == semiring->add->op->ztype) ;
66 GB_Opcode opcode = mult->opcode ;
67 // GB_reduce_to_vector does not use GB_AxB_rowscale:
68 ASSERT (!(mult->function == NULL &&
69 (opcode == GB_FIRST_opcode || opcode == GB_SECOND_opcode))) ;
70
71 //--------------------------------------------------------------------------
72 // copy the pattern of B into C
73 //--------------------------------------------------------------------------
74
75 // allocate C->x but do not initialize it
76 info = GB_dup2 (&C, B, false, mult->ztype, Context) ; // static header
77 if (info != GrB_SUCCESS)
78 {
79 // out of memory
80 return (info) ;
81 }
82
83 //--------------------------------------------------------------------------
84 // apply a positional operator: convert C=D*B to C=op(B)
85 //--------------------------------------------------------------------------
86
87 if (GB_OPCODE_IS_POSITIONAL (opcode))
88 {
89 if (flipxy)
90 {
91 // the multiplicatve operator is fmult(y,x), so flip the opcode
92 bool handled ;
93 opcode = GB_binop_flip (opcode, &handled) ; // for positional ops
94 ASSERT (handled) ; // all positional ops can be flipped
95 }
96 // determine unary operator to compute C=D*B
97 GrB_UnaryOp op1 = NULL ;
98 if (mult->ztype == GrB_INT64)
99 {
100 switch (opcode)
101 {
102 // first_op(D,B) becomes position_i(B)
103 case GB_FIRSTI_opcode :
104 case GB_FIRSTJ_opcode : op1 = GxB_POSITIONI_INT64 ;
105 break ;
106 case GB_FIRSTI1_opcode :
107 case GB_FIRSTJ1_opcode : op1 = GxB_POSITIONI1_INT64 ;
108 break ;
109 // second_op(D,B) becomes position_op(B)
110 case GB_SECONDI_opcode : op1 = GxB_POSITIONI_INT64 ;
111 break ;
112 case GB_SECONDJ_opcode : op1 = GxB_POSITIONJ_INT64 ;
113 break ;
114 case GB_SECONDI1_opcode : op1 = GxB_POSITIONI1_INT64 ;
115 break ;
116 case GB_SECONDJ1_opcode : op1 = GxB_POSITIONJ1_INT64 ;
117 break ;
118 default: ;
119 }
120 }
121 else
122 {
123 switch (opcode)
124 {
125 // first_op(D,B) becomes position_i(B)
126 case GB_FIRSTI_opcode :
127 case GB_FIRSTJ_opcode : op1 = GxB_POSITIONI_INT32 ;
128 break ;
129 case GB_FIRSTI1_opcode :
130 case GB_FIRSTJ1_opcode : op1 = GxB_POSITIONI1_INT32 ;
131 break ;
132 // second_op(D,B) becomes position_op(B)
133 case GB_SECONDI_opcode : op1 = GxB_POSITIONI_INT32 ;
134 break ;
135 case GB_SECONDJ_opcode : op1 = GxB_POSITIONJ_INT32 ;
136 break ;
137 case GB_SECONDI1_opcode : op1 = GxB_POSITIONI1_INT32 ;
138 break ;
139 case GB_SECONDJ1_opcode : op1 = GxB_POSITIONJ1_INT32 ;
140 break ;
141 default: ;
142 }
143 }
144 GB_OK (GB_apply_op ((GB_void *) (C->x), op1, // positional unary op
145 NULL, NULL, false, B, Context)) ;
146 ASSERT_MATRIX_OK (C, "rowscale positional: C = D*B output", GB0) ;
147 return (GrB_SUCCESS) ;
148 }
149
150 //--------------------------------------------------------------------------
151 // determine the number of threads to use
152 //--------------------------------------------------------------------------
153
154 GB_GET_NTHREADS_MAX (nthreads_max, chunk, Context) ;
155 int nthreads = GB_nthreads (GB_NNZ_HELD (B) + B->nvec, chunk, nthreads_max);
156
157 //--------------------------------------------------------------------------
158 // determine if the values are accessed
159 //--------------------------------------------------------------------------
160
161 bool op_is_first = (opcode == GB_FIRST_opcode) ;
162 bool op_is_second = (opcode == GB_SECOND_opcode) ;
163 bool op_is_pair = (opcode == GB_PAIR_opcode) ;
164 bool D_is_pattern = false ;
165 bool B_is_pattern = false ;
166
167 if (flipxy)
168 {
169 // z = fmult (b,a) will be computed
170 D_is_pattern = op_is_first || op_is_pair ;
171 B_is_pattern = op_is_second || op_is_pair ;
172 ASSERT (GB_IMPLIES (!D_is_pattern,
173 GB_Type_compatible (D->type, mult->ytype))) ;
174 ASSERT (GB_IMPLIES (!B_is_pattern,
175 GB_Type_compatible (B->type, mult->xtype))) ;
176 }
177 else
178 {
179 // z = fmult (a,b) will be computed
180 D_is_pattern = op_is_second || op_is_pair ;
181 B_is_pattern = op_is_first || op_is_pair ;
182 ASSERT (GB_IMPLIES (!D_is_pattern,
183 GB_Type_compatible (D->type, mult->xtype))) ;
184 ASSERT (GB_IMPLIES (!B_is_pattern,
185 GB_Type_compatible (B->type, mult->ytype))) ;
186 }
187
188 //--------------------------------------------------------------------------
189 // C = D*B, row scale, via built-in binary operators
190 //--------------------------------------------------------------------------
191
192 bool done = false ;
193
194 #ifndef GBCOMPACT
195
196 //----------------------------------------------------------------------
197 // define the worker for the switch factory
198 //----------------------------------------------------------------------
199
200 #define GB_DxB(mult,xname) GB (_DxB_ ## mult ## xname)
201
202 #define GB_BINOP_WORKER(mult,xname) \
203 { \
204 info = GB_DxB(mult,xname) (C, D, D_is_pattern, B, B_is_pattern, \
205 nthreads) ; \
206 done = (info != GrB_NO_VALUE) ; \
207 } \
208 break ;
209
210 //----------------------------------------------------------------------
211 // launch the switch factory
212 //----------------------------------------------------------------------
213
214 GB_Type_code xcode, ycode, zcode ;
215 if (GB_binop_builtin (D->type, D_is_pattern, B->type, B_is_pattern,
216 mult, flipxy, &opcode, &xcode, &ycode, &zcode))
217 {
218 // C=D*B, rowscale with built-in operator
219 #define GB_BINOP_IS_SEMIRING_MULTIPLIER
220 #include "GB_binop_factory.c"
221 #undef GB_BINOP_IS_SEMIRING_MULTIPLIER
222 }
223
224 #endif
225
226 //--------------------------------------------------------------------------
227 // C = D*B, row scale, with typecasting or user-defined operator
228 //--------------------------------------------------------------------------
229
230 if (!done)
231 {
232
233 //----------------------------------------------------------------------
234 // get operators, functions, workspace, contents of D, B, and C
235 //----------------------------------------------------------------------
236
237 GB_BURBLE_MATRIX (C, "(generic C=D*B rowscale) ") ;
238
239 GxB_binary_function fmult = mult->function ;
240
241 size_t csize = C->type->size ;
242 size_t dsize = D_is_pattern ? 0 : D->type->size ;
243 size_t bsize = B_is_pattern ? 0 : B->type->size ;
244
245 size_t xsize = mult->xtype->size ;
246 size_t ysize = mult->ytype->size ;
247
248 // scalar workspace: because of typecasting, the x/y types need not
249 // be the same as the size of the D and B types.
250 // flipxy false: dii = (xtype) D(i,i) and bij = (ytype) B(i,j)
251 // flipxy true: dii = (ytype) D(i,i) and bij = (xtype) B(i,j)
252 size_t dii_size = flipxy ? ysize : xsize ;
253 size_t bij_size = flipxy ? xsize : ysize ;
254
255 GB_void *restrict Cx = (GB_void *) C->x ;
256
257 GB_cast_function cast_D, cast_B ;
258 if (flipxy)
259 {
260 // D is typecasted to y, and B is typecasted to x
261 cast_D = D_is_pattern ? NULL :
262 GB_cast_factory (mult->ytype->code, D->type->code) ;
263 cast_B = B_is_pattern ? NULL :
264 GB_cast_factory (mult->xtype->code, B->type->code) ;
265 }
266 else
267 {
268 // D is typecasted to x, and B is typecasted to y
269 cast_D = D_is_pattern ? NULL :
270 GB_cast_factory (mult->xtype->code, D->type->code) ;
271 cast_B = B_is_pattern ? NULL :
272 GB_cast_factory (mult->ytype->code, B->type->code) ;
273 }
274
275 //----------------------------------------------------------------------
276 // C = D*B via function pointers, and typecasting
277 //----------------------------------------------------------------------
278
279 // dii = D(i,i), located in Dx [i]
280 #define GB_GETA(dii,Dx,i) \
281 GB_void dii [GB_VLA(dii_size)] ; \
282 if (!D_is_pattern) cast_D (dii, Dx +((i)*dsize), dsize) ;
283
284 // bij = B(i,j), located in Bx [pB]
285 #define GB_GETB(bij,Bx,pB) \
286 GB_void bij [GB_VLA(bij_size)] ; \
287 if (!B_is_pattern) cast_B (bij, Bx +((pB)*bsize), bsize) ;
288
289 // address of Cx [p]
290 #define GB_CX(p) Cx +((p)*csize)
291
292 #define GB_ATYPE GB_void
293 #define GB_BTYPE GB_void
294 #define GB_CTYPE GB_void
295
296 // no vectorization
297 #define GB_PRAGMA_SIMD_VECTORIZE ;
298
299 if (flipxy)
300 {
301 #define GB_BINOP(z,x,y,i,j) fmult (z,y,x)
302 #include "GB_AxB_rowscale_meta.c"
303 #undef GB_BINOP
304 }
305 else
306 {
307 #define GB_BINOP(z,x,y,i,j) fmult (z,x,y)
308 #include "GB_AxB_rowscale_meta.c"
309 #undef GB_BINOP
310 }
311 }
312
313 //--------------------------------------------------------------------------
314 // return result
315 //--------------------------------------------------------------------------
316
317 ASSERT_MATRIX_OK (C, "rowscale: C = D*B output", GB0) ;
318 return (GrB_SUCCESS) ;
319 }
320
321