1 //------------------------------------------------------------------------------
2 // GB_mex_rdiv: compute C=A*B with the rdiv 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-rdiv-fp64 where plus is the
12 // built-in GrB_PLUS_FP64 operator, and rdiv is z=y/x in double.
13
14 #include "GB_mex.h"
15
16 #define USAGE "C = GB_mex_rdiv (A, B, axb_method, cprint)"
17
18 #define FREE_ALL \
19 { \
20 GrB_Matrix_free_(&A) ; \
21 GrB_Matrix_free_(&B) ; \
22 GrB_Matrix_free_(&C) ; \
23 GrB_BinaryOp_free_(&My_rdiv) ; \
24 GrB_Semiring_free_(&My_plus_rdiv) ; \
25 GB_mx_put_global (true) ; \
26 }
27
28 //------------------------------------------------------------------------------
29
30 GrB_Info info ;
31 bool malloc_debug = false ;
32 bool ignore = false, ignore1 = false, ignore2 = false ;
33 bool cprint = false ;
34 GrB_Matrix A = NULL, B = NULL, C = NULL, MT = NULL ;
35 int64_t anrows = 0 ;
36 int64_t ancols = 0 ;
37 int64_t bnrows = 0 ;
38 int64_t bncols = 0 ;
39 GrB_Desc_Value AxB_method = GxB_DEFAULT ;
40 struct GB_Matrix_opaque MT_header, C_header ;
41
42 GrB_Info axb (GB_Context Context, bool cprint) ;
43
44 GrB_Semiring My_plus_rdiv = NULL ;
45 GrB_BinaryOp My_rdiv = NULL ;
46
47 void my_rdiv (double *z, const double *x, const double *y) ;
48
my_rdiv(double * z,const double * x,const double * y)49 void my_rdiv (double *z, const double *x, const double *y)
50 {
51 (*z) = (*y) / (*x) ;
52 }
53
54 //------------------------------------------------------------------------------
55
axb(GB_Context Context,bool cprint)56 GrB_Info axb (GB_Context Context, bool cprint)
57 {
58 // create the rdiv operator
59 info = GrB_BinaryOp_new (&My_rdiv, my_rdiv, GrB_FP64, GrB_FP64, GrB_FP64) ;
60 if (info != GrB_SUCCESS) return (info) ;
61 GrB_BinaryOp_wait_(&My_rdiv) ;
62 if (info != GrB_SUCCESS) return (info) ;
63 info = GrB_Semiring_new (&My_plus_rdiv, GxB_PLUS_FP64_MONOID, My_rdiv) ;
64 if (info != GrB_SUCCESS)
65 {
66 GrB_BinaryOp_free_(&My_rdiv) ;
67 return (info) ;
68 }
69
70 MT = GB_clear_static_header (&MT_header) ;
71 C = GB_clear_static_header (&C_header) ;
72
73 // C = A*B
74 info = GB_AxB_meta (C, NULL, // C cannot be computed in place
75 false, // C_replace
76 true, // CSC
77 MT, // no MT returned
78 &ignore1, // M_transposed will be false
79 NULL, // no Mask
80 false, // mask not complemented
81 false, // mask not structural
82 NULL, // no accum
83 A, B,
84 My_plus_rdiv,
85 false, // no A transpose
86 false, // no B transpose
87 false, // no flipxy
88 &ignore, // mask_applied
89 &ignore2, // done_in_place
90 AxB_method,
91 true, // do the sort
92 Context) ;
93
94 if (C != NULL)
95 {
96 if (cprint) GxB_Matrix_fprint_(C, GxB_COMPLETE, NULL) ;
97 }
98
99 GrB_BinaryOp_free_(&My_rdiv) ;
100 GrB_Semiring_free_(&My_plus_rdiv) ;
101
102 return (info) ;
103 }
104
105 //------------------------------------------------------------------------------
106
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])107 void mexFunction
108 (
109 int nargout,
110 mxArray *pargout [ ],
111 int nargin,
112 const mxArray *pargin [ ]
113 )
114 {
115
116 info = GrB_SUCCESS ;
117 malloc_debug = GB_mx_get_global (true) ;
118 ignore = false ;
119 ignore1 = false ;
120 ignore2 = false ;
121 A = NULL ;
122 B = NULL ;
123 C = NULL ;
124
125 My_rdiv = NULL ;
126 My_plus_rdiv = NULL ;
127
128 GB_CONTEXT (USAGE) ;
129
130 // check inputs
131 if (nargout > 1 || nargin < 2 || nargin > 4)
132 {
133 mexErrMsgTxt ("Usage: " USAGE) ;
134 }
135
136 #define GET_DEEP_COPY ;
137 #define FREE_DEEP_COPY ;
138
139 // get A and B
140 A = GB_mx_mxArray_to_Matrix (pargin [0], "A", false, true) ;
141 B = GB_mx_mxArray_to_Matrix (pargin [1], "B", false, true) ;
142 if (A == NULL || B == NULL)
143 {
144 FREE_ALL ;
145 mexErrMsgTxt ("failed") ;
146 }
147
148 if (!A->is_csc || !B->is_csc)
149 {
150 mexErrMsgTxt ("A and B must be in CSC format") ;
151 }
152
153 // get the axb_method
154 // 0 or not present: default
155 // 1001: Gustavson
156 // 1003: dot
157 // 1004: hash
158 // 1005: saxpy
159 GET_SCALAR (2, GrB_Desc_Value, AxB_method, GxB_DEFAULT) ;
160
161 // get the cprint flag
162 GET_SCALAR (3, bool, cprint, false) ;
163
164 if (! ((AxB_method == GxB_DEFAULT) ||
165 (AxB_method == GxB_AxB_GUSTAVSON) ||
166 (AxB_method == GxB_AxB_HASH) ||
167 (AxB_method == GxB_AxB_SAXPY) ||
168 (AxB_method == GxB_AxB_DOT)))
169 {
170 mexErrMsgTxt ("unknown method") ;
171 }
172
173 // determine the dimensions
174 anrows = GB_NROWS (A) ;
175 ancols = GB_NCOLS (A) ;
176 bnrows = GB_NROWS (B) ;
177 bncols = GB_NCOLS (B) ;
178 if (ancols != bnrows)
179 {
180 FREE_ALL ;
181 mexErrMsgTxt ("invalid dimensions") ;
182 }
183
184 METHOD (axb (Context, cprint)) ;
185
186 // return C to MATLAB
187 pargout [0] = GB_mx_Matrix_to_mxArray (&C, "C AxB result", false) ;
188
189 FREE_ALL ;
190 }
191
192