1 //------------------------------------------------------------------------------
2 // GB_apply: apply a unary operator; optionally transpose a matrix
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, op(A)) or accum (C, op(A)')
11
12 // GB_apply does the work for GrB_*_apply, including the binary op variants.
13
14 #include "GB_apply.h"
15 #include "GB_transpose.h"
16 #include "GB_accum_mask.h"
17
18 #define GB_FREE_ALL ;
19
GB_apply(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_UnaryOp op1_in,const GrB_BinaryOp op2_in,const GxB_Scalar scalar,bool binop_bind1st,const GrB_Matrix A,bool A_transpose,GB_Context Context)20 GrB_Info GB_apply // C<M> = accum (C, op(A)) or op(A')
21 (
22 GrB_Matrix C, // input/output matrix for results
23 const bool C_replace, // C descriptor
24 const GrB_Matrix M, // optional mask for C, unused if NULL
25 const bool Mask_comp, // M descriptor
26 const bool Mask_struct, // if true, use the only structure of M
27 const GrB_BinaryOp accum, // optional accum for Z=accum(C,T)
28 const GrB_UnaryOp op1_in, // unary operator to apply
29 const GrB_BinaryOp op2_in, // binary operator to apply
30 const GxB_Scalar scalar, // scalar to bind to binary operator
31 bool binop_bind1st, // if true, binop(x,A) else binop(A,y)
32 const GrB_Matrix A, // first input: matrix A
33 bool A_transpose, // A matrix descriptor
34 GB_Context Context
35 )
36 {
37
38 //--------------------------------------------------------------------------
39 // check inputs
40 //--------------------------------------------------------------------------
41
42 // C may be aliased with M and/or A
43
44 struct GB_Matrix_opaque T_header ;
45 GrB_Matrix T = GB_clear_static_header (&T_header) ;
46
47 GB_RETURN_IF_FAULTY_OR_POSITIONAL (accum) ;
48 ASSERT_MATRIX_OK (C, "C input for GB_apply", GB0) ;
49 ASSERT_MATRIX_OK_OR_NULL (M, "M for GB_apply", GB0) ;
50 ASSERT_BINARYOP_OK_OR_NULL (accum, "accum for GB_apply", GB0) ;
51 ASSERT_MATRIX_OK (A, "A input for GB_apply", GB0) ;
52
53 GrB_UnaryOp op1 = op1_in ;
54 GrB_BinaryOp op2 = op2_in ;
55 GB_Opcode opcode ;
56 GrB_Type T_type ;
57 if (op1 != NULL)
58 {
59 // apply a unary operator
60 GB_RETURN_IF_FAULTY (op1) ;
61 ASSERT_UNARYOP_OK (op1, "op1 for GB_apply", GB0) ;
62 T_type = op1->ztype ;
63 opcode = op1->opcode ;
64 if (!GB_OPCODE_IS_POSITIONAL (opcode))
65 {
66 // A must also be compatible with op1->xtype
67 if (!GB_Type_compatible (A->type, op1->xtype))
68 {
69 GB_ERROR (GrB_DOMAIN_MISMATCH,
70 "Incompatible type for z=%s(x):\n"
71 "input A of type [%s]\n"
72 "cannot be typecast to x input of type [%s]",
73 op1->name, A->type->name, op1->xtype->name) ;
74 }
75 }
76 }
77 else if (op2 != NULL)
78 {
79 // apply a binary operator, with one input bound to a scalar
80 GB_RETURN_IF_FAULTY (op2) ;
81 ASSERT_BINARYOP_OK (op2, "op2 for GB_apply", GB0) ;
82 ASSERT_SCALAR_OK (scalar, "scalar for GB_apply", GB0) ;
83 T_type = op2->ztype ;
84 opcode = op2->opcode ;
85 if (!GB_OPCODE_IS_POSITIONAL (opcode))
86 {
87 bool op_is_first = opcode == GB_FIRST_opcode ;
88 bool op_is_second = opcode == GB_SECOND_opcode ;
89 bool op_is_pair = opcode == GB_PAIR_opcode ;
90 if (binop_bind1st)
91 {
92 // C = op (scalar,A)
93 // A must be compatible with op2->ytype
94 if (!(op_is_first || op_is_pair ||
95 GB_Type_compatible (A->type, op2->ytype)))
96 {
97 GB_ERROR (GrB_DOMAIN_MISMATCH,
98 "Incompatible type for z=%s(x,y):\n"
99 "input A of type [%s]\n"
100 "cannot be typecast to y input of type [%s]",
101 op2->name, A->type->name, op2->ytype->name) ;
102 }
103 // scalar must be compatible with op2->xtype
104 if (!(op_is_second || op_is_pair ||
105 GB_Type_compatible (scalar->type, op2->xtype)))
106 {
107 GB_ERROR (GrB_DOMAIN_MISMATCH,
108 "Incompatible type for z=%s(x,y):\n"
109 "input scalar of type [%s]\n"
110 "cannot be typecast to x input of type [%s]",
111 op2->name, scalar->type->name, op2->xtype->name) ;
112 }
113 }
114 else
115 {
116 // C = op (A,scalar)
117 // A must be compatible with op2->xtype
118 if (!(op_is_first || op_is_pair ||
119 GB_Type_compatible (A->type, op2->xtype)))
120 {
121 GB_ERROR (GrB_DOMAIN_MISMATCH,
122 "Incompatible type for z=%s(x,y):\n"
123 "input scalar of type [%s]\n"
124 "cannot be typecast to x input of type [%s]",
125 op2->name, A->type->name, op2->xtype->name) ;
126 }
127 // scalar must be compatible with op2->ytype
128 if (!(op_is_second || op_is_pair
129 || GB_Type_compatible (scalar->type, op2->ytype)))
130 {
131 GB_ERROR (GrB_DOMAIN_MISMATCH,
132 "Incompatible type for z=%s(x,y):\n"
133 "input A of type [%s]\n"
134 "cannot be typecast to y input of type [%s]",
135 op2->name, scalar->type->name, op2->ytype->name) ;
136 }
137 }
138 }
139 }
140 else
141 {
142 GB_ERROR (GrB_NULL_POINTER,
143 "Required argument is null: [%s]", "op") ;
144 }
145
146 // check domains and dimensions for C<M> = accum (C,T)
147 GrB_Info info ;
148 GB_OK (GB_compatible (C->type, C, M, Mask_struct, accum, T_type, Context)) ;
149
150 // check the dimensions
151 int64_t tnrows = (A_transpose) ? GB_NCOLS (A) : GB_NROWS (A) ;
152 int64_t tncols = (A_transpose) ? GB_NROWS (A) : GB_NCOLS (A) ;
153 if (GB_NROWS (C) != tnrows || GB_NCOLS (C) != tncols)
154 {
155 GB_ERROR (GrB_DIMENSION_MISMATCH,
156 "Dimensions not compatible:\n"
157 "output is " GBd "-by-" GBd "\n"
158 "input is " GBd "-by-" GBd "%s",
159 GB_NROWS (C), GB_NCOLS (C),
160 tnrows, tncols, A_transpose ? " (transposed)" : "") ;
161 }
162
163 // quick return if an empty mask is complemented
164 GB_RETURN_IF_QUICK_MASK (C, C_replace, M, Mask_comp, Mask_struct) ;
165
166 // delete any lingering zombies and assemble any pending tuples
167 GB_MATRIX_WAIT_IF_PENDING_OR_ZOMBIES (A) ; // A can be jumbled
168 GB_MATRIX_WAIT (scalar) ;
169
170 if (op2 != NULL && GB_NNZ (scalar) != 1)
171 {
172 // the scalar entry must be present
173 GB_ERROR (GrB_INVALID_VALUE, "%s", "Scalar must contain an entry") ;
174 }
175
176 //--------------------------------------------------------------------------
177 // rename first, second, any, and pair operators
178 //--------------------------------------------------------------------------
179
180 if (op2 != NULL)
181 {
182 // first(A,x), second(y,A), and any(...) become identity(A)
183 if ((opcode == GB_ANY_opcode) ||
184 (opcode == GB_FIRST_opcode && !binop_bind1st) ||
185 (opcode == GB_SECOND_opcode && binop_bind1st))
186 {
187 switch (op2->xtype->code)
188 {
189 default :
190 case GB_BOOL_code : op1 = GrB_IDENTITY_BOOL ; break ;
191 case GB_INT8_code : op1 = GrB_IDENTITY_INT8 ; break ;
192 case GB_INT16_code : op1 = GrB_IDENTITY_INT16 ; break ;
193 case GB_INT32_code : op1 = GrB_IDENTITY_INT32 ; break ;
194 case GB_INT64_code : op1 = GrB_IDENTITY_INT64 ; break ;
195 case GB_UINT8_code : op1 = GrB_IDENTITY_UINT8 ; break ;
196 case GB_UINT16_code : op1 = GrB_IDENTITY_UINT16 ; break ;
197 case GB_UINT32_code : op1 = GrB_IDENTITY_UINT32 ; break ;
198 case GB_UINT64_code : op1 = GrB_IDENTITY_UINT64 ; break ;
199 case GB_FP32_code : op1 = GrB_IDENTITY_FP32 ; break ;
200 case GB_FP64_code : op1 = GrB_IDENTITY_FP64 ; break ;
201 case GB_FC32_code : op1 = GxB_IDENTITY_FC32 ; break ;
202 case GB_FC64_code : op1 = GxB_IDENTITY_FC64 ; break ;
203 }
204 op2 = NULL ;
205 }
206 else if (opcode == GB_PAIR_opcode)
207 {
208 // pair (...) becomes one(A)
209 switch (op2->xtype->code)
210 {
211 default :
212 case GB_BOOL_code : op1 = GxB_ONE_BOOL ; break ;
213 case GB_INT8_code : op1 = GxB_ONE_INT8 ; break ;
214 case GB_INT16_code : op1 = GxB_ONE_INT16 ; break ;
215 case GB_INT32_code : op1 = GxB_ONE_INT32 ; break ;
216 case GB_INT64_code : op1 = GxB_ONE_INT64 ; break ;
217 case GB_UINT8_code : op1 = GxB_ONE_UINT8 ; break ;
218 case GB_UINT16_code : op1 = GxB_ONE_UINT16 ; break ;
219 case GB_UINT32_code : op1 = GxB_ONE_UINT32 ; break ;
220 case GB_UINT64_code : op1 = GxB_ONE_UINT64 ; break ;
221 case GB_FP32_code : op1 = GxB_ONE_FP32 ; break ;
222 case GB_FP64_code : op1 = GxB_ONE_FP64 ; break ;
223 case GB_FC32_code : op1 = GxB_ONE_FC32 ; break ;
224 case GB_FC64_code : op1 = GxB_ONE_FC64 ; break ;
225 }
226 op2 = NULL ;
227 }
228 }
229
230 //--------------------------------------------------------------------------
231 // T = op(A) or op(A')
232 //--------------------------------------------------------------------------
233
234 bool T_is_csc = C->is_csc ;
235 if (T_is_csc != A->is_csc)
236 {
237 // Flip the sense of A_transpose
238 A_transpose = !A_transpose ;
239 }
240
241 if (!T_is_csc)
242 {
243 // positional ops must be flipped, with i and j swapped
244 if (op1 != NULL)
245 {
246 op1 = GB_positional_unop_ijflip (op1) ;
247 opcode = op1->opcode ;
248 }
249 else if (op2 != NULL)
250 {
251 op2 = GB_positional_binop_ijflip (op2) ;
252 opcode = op2->opcode ;
253 }
254 }
255
256 if (A_transpose)
257 {
258 // T = op (A'), typecasting to op*->ztype
259 // transpose: typecast, apply an op, not in-place.
260 GBURBLE ("(transpose-op) ") ;
261 info = GB_transpose (&T, T_type, T_is_csc, A, // T static
262 op1, op2, scalar, binop_bind1st, Context) ;
263 ASSERT (GB_JUMBLED_OK (T)) ;
264 // A positional op is applied to C after the transpose is computed,
265 // using the T_is_csc format. The ijflip is handled
266 // above.
267 }
268 else if (M == NULL && accum == NULL && (C == A) && C->type == T_type)
269 {
270 GBURBLE ("(inplace-op) ") ;
271 // C = op (C), operating on the values in-place, with no typecasting
272 // of the output of the operator with the matrix C.
273 // No work to do if the op is identity.
274 // FUTURE::: also handle C += op(C), with accum.
275 if (opcode != GB_IDENTITY_opcode)
276 {
277 // the output Cx is aliased with C->x in GB_apply_op.
278 GB_void *Cx = (GB_void *) C->x ;
279 info = GB_apply_op (Cx, op1, op2, // op1 != identity
280 scalar, binop_bind1st, C, Context) ;
281 }
282 return (info) ;
283 }
284 else
285 {
286 // T = op (A), pattern is a shallow copy of A, type is op*->ztype.
287 GBURBLE ("(shallow-op) ") ;
288 info = GB_shallow_op (T, T_is_csc,
289 op1, op2, scalar, binop_bind1st, A, Context) ;
290 }
291
292 if (info != GrB_SUCCESS)
293 {
294 GB_phbix_free (T) ;
295 return (info) ;
296 }
297
298 ASSERT (T->is_csc == C->is_csc) ;
299
300 //--------------------------------------------------------------------------
301 // C<M> = accum (C,T): accumulate the results into C via the M
302 //--------------------------------------------------------------------------
303
304 return (GB_accum_mask (C, M, NULL, accum, &T, C_replace, Mask_comp,
305 Mask_struct, Context)) ;
306 }
307
308