1 //------------------------------------------------------------------------------
2 // GB_mex_rdiv2: compute C=A*B with the rdiv2 operator
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 // This is for testing only. See GrB_mxm instead. Returns a plain MATLAB
11 // matrix, in double. The semiring is plus-rdiv2 where plus is the
12 // built-in GrB_PLUS_FP64 operator, and rdiv2 is z=y/x with y float and x
13 // double. The input matrix B is typecasted here to GrB_FP32.
14
15 #include "GB_mex.h"
16
17 #define USAGE "C = GB_mex_rdiv2 (A, B, atrans, btrans, axb_method, flipxy, C_scalar)"
18
19 #define FREE_ALL \
20 { \
21 GrB_Matrix_free_(&A) ; \
22 GrB_Matrix_free_(&B) ; \
23 GrB_Matrix_free_(&B64) ; \
24 GrB_Matrix_free_(&C) ; \
25 GrB_Matrix_free_(&T) ; \
26 GrB_BinaryOp_free_(&My_rdiv2) ; \
27 GrB_Semiring_free_(&My_plus_rdiv2) ; \
28 GB_mx_put_global (true) ; \
29 }
30
31 //------------------------------------------------------------------------------
32
33 GrB_Info info ;
34 bool malloc_debug = false ;
35 bool ignore = false, ignore1 = false, ignore2 = false ;
36 bool atranspose = false ;
37 bool btranspose = false ;
38 GrB_Matrix A = NULL, B = NULL, B64 = NULL, C = NULL, T = NULL, MT = NULL ;
39 int64_t anrows = 0 ;
40 int64_t ancols = 0 ;
41 int64_t bnrows = 0 ;
42 int64_t bncols = 0 ;
43 GrB_Desc_Value AxB_method = GxB_DEFAULT ;
44 bool flipxy = false ;
45 bool done_in_place = false ;
46 double C_scalar = 0 ;
47 struct GB_Matrix_opaque MT_header, T_header ;
48
49 GrB_Info axb (GB_Context Context) ;
50
51 GrB_Semiring My_plus_rdiv2 = NULL ;
52 GrB_BinaryOp My_rdiv2 = NULL ;
53
54 void my_rdiv2 (double *z, const double *x, const float *y) ;
55
my_rdiv2(double * z,const double * x,const float * y)56 void my_rdiv2 (double *z, const double *x, const float *y)
57 {
58 (*z) = ((double) (*y)) / (*x) ;
59 }
60
61 //------------------------------------------------------------------------------
62
axb(GB_Context Context)63 GrB_Info axb (GB_Context Context)
64 {
65 // create the rdiv2 operator
66 info = GrB_BinaryOp_new (&My_rdiv2, my_rdiv2, GrB_FP64, GrB_FP64, GrB_FP32);
67 GrB_BinaryOp_wait_(&My_rdiv2) ;
68 if (info != GrB_SUCCESS) return (info) ;
69 info = GrB_Semiring_new (&My_plus_rdiv2, GxB_PLUS_FP64_MONOID, My_rdiv2) ;
70 if (info != GrB_SUCCESS)
71 {
72 GrB_BinaryOp_free_(&My_rdiv2) ;
73 return (info) ;
74 }
75
76 bool do_in_place = (C_scalar != 0) ;
77 C = NULL ;
78
79 if (do_in_place)
80 {
81 // construct the result matrix and fill it with the scalar
82 GrB_Index cnrows = anrows ;
83 GrB_Index cncols = bncols ;
84 info = GrB_Matrix_new (&C, GrB_FP64, cnrows, cncols) ;
85 if (info != GrB_SUCCESS)
86 {
87 GrB_BinaryOp_free_(&My_rdiv2) ;
88 GrB_Semiring_free_(&My_plus_rdiv2) ;
89 return (info) ;
90 }
91 info = GrB_Matrix_assign_FP64_(C, NULL, NULL, C_scalar,
92 GrB_ALL, cnrows, GrB_ALL, cncols, NULL) ;
93 if (info != GrB_SUCCESS)
94 {
95 GrB_BinaryOp_free_(&My_rdiv2) ;
96 GrB_Semiring_free_(&My_plus_rdiv2) ;
97 GrB_Matrix_free_(&C) ;
98 return (info) ;
99 }
100 }
101
102 MT = GB_clear_static_header (&MT_header) ;
103 T = GB_clear_static_header (&T_header) ;
104
105 // C = A*B or C += A*B
106 info = GB_AxB_meta (T, C, // can be done in place if C != NULL
107 false, // C_replace
108 true, // CSC
109 MT, // no MT returned
110 &ignore1, // M_transposed will be false
111 NULL, // no Mask
112 false, // mask not complemented
113 false, // mask not structural
114 (do_in_place) ? GrB_PLUS_FP64 : NULL, // accum
115 A, B,
116 My_plus_rdiv2,
117 atranspose,
118 btranspose,
119 flipxy,
120 &ignore, // mask_applied
121 &done_in_place,
122 AxB_method,
123 true, // do the sort
124 Context) ;
125
126 if (info == GrB_SUCCESS)
127 {
128 if (done_in_place != do_in_place)
129 {
130 mexErrMsgTxt ("failure: not in place as expected\n") ;
131 }
132 if (!done_in_place)
133 {
134 GrB_Matrix_free_(&C) ;
135 info = GrB_Matrix_dup (&C, T) ;
136 }
137 }
138
139 if (info != GrB_SUCCESS)
140 {
141 GrB_Matrix_free_(&C) ;
142 }
143
144 GrB_Matrix_free_(&T) ;
145
146 GrB_BinaryOp_free_(&My_rdiv2) ;
147 GrB_Semiring_free_(&My_plus_rdiv2) ;
148
149 return (info) ;
150 }
151
152 //------------------------------------------------------------------------------
153
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])154 void mexFunction
155 (
156 int nargout,
157 mxArray *pargout [ ],
158 int nargin,
159 const mxArray *pargin [ ]
160 )
161 {
162
163 info = GrB_SUCCESS ;
164 malloc_debug = GB_mx_get_global (true) ;
165 ignore = false ;
166 ignore1 = false ;
167 ignore2 = false ;
168 A = NULL ;
169 B = NULL ;
170 B64 = NULL ;
171 C = NULL ;
172
173 My_rdiv2 = NULL ;
174 My_plus_rdiv2 = NULL ;
175
176 GB_CONTEXT (USAGE) ;
177
178 // check inputs
179 if (nargout > 1 || nargin < 2 || nargin > 7)
180 {
181 mexErrMsgTxt ("Usage: " USAGE) ;
182 }
183
184 #define GET_DEEP_COPY ;
185 #define FREE_DEEP_COPY ;
186
187 // get A and B64
188 A = GB_mx_mxArray_to_Matrix (pargin [0], "A", false, true) ;
189 B64 = GB_mx_mxArray_to_Matrix (pargin [1], "B", false, true) ;
190 if (A == NULL || B64 == NULL)
191 {
192 FREE_ALL ;
193 mexErrMsgTxt ("failed") ;
194 }
195
196 if (!A->is_csc || !B64->is_csc)
197 {
198 mexErrMsgTxt ("A and B must be in CSC format") ;
199 }
200
201 // get the atranspose option
202 GET_SCALAR (2, bool, atranspose, false) ;
203
204 // get the btranspose option
205 GET_SCALAR (3, bool, btranspose, false) ;
206
207 // get the axb_method
208 // 0 or not present: default
209 // 1001: Gustavson
210 // 1003: dot
211 // 1004: hash
212 // 1005: saxpy
213 GET_SCALAR (4, GrB_Desc_Value, AxB_method, GxB_DEFAULT) ;
214
215 if (! ((AxB_method == GxB_DEFAULT) ||
216 (AxB_method == GxB_AxB_GUSTAVSON) ||
217 (AxB_method == GxB_AxB_HASH) ||
218 (AxB_method == GxB_AxB_SAXPY) ||
219 (AxB_method == GxB_AxB_DOT)))
220 {
221 mexErrMsgTxt ("unknown method") ;
222 }
223
224 // get the flipxy option
225 GET_SCALAR (5, bool, flipxy, false) ;
226
227 // get the C_scalar
228 GET_SCALAR (6, double, C_scalar, 0) ;
229
230 // determine the dimensions
231 anrows = (atranspose) ? GB_NCOLS (A) : GB_NROWS (A) ;
232 ancols = (atranspose) ? GB_NROWS (A) : GB_NCOLS (A) ;
233 bnrows = (btranspose) ? GB_NCOLS (B64) : GB_NROWS (B64) ;
234 bncols = (btranspose) ? GB_NROWS (B64) : GB_NCOLS (B64) ;
235 if (ancols != bnrows)
236 {
237 FREE_ALL ;
238 mexErrMsgTxt ("invalid dimensions") ;
239 }
240
241 if (atranspose && btranspose && C_scalar != 0)
242 {
243 C_scalar = 0 ;
244 }
245
246 // convert B64 (double) to B (float)
247 GrB_Matrix_new (&B, GrB_FP32, bnrows, bncols) ;
248 GrB_Matrix_assign_(B, NULL, NULL, B64, GrB_ALL, 0, GrB_ALL, 0, NULL) ;
249
250 // B must be completed
251 GrB_Matrix_wait (&B) ;
252
253 METHOD (axb (Context)) ;
254
255 // return C to MATLAB
256 pargout [0] = GB_mx_Matrix_to_mxArray (&C, "C AxB result", false) ;
257
258 FREE_ALL ;
259 }
260
261