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
|
//------------------------------------------------------------------------------
// GB_mex_rdiv: compute C=A*B with the rdiv operator
//------------------------------------------------------------------------------
// SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2022, All Rights Reserved.
// SPDX-License-Identifier: Apache-2.0
//------------------------------------------------------------------------------
// This is for testing only. See GrB_mxm instead. Returns a plain built-in
// matrix, in double. The semiring is plus-rdiv-fp64 where plus is the
// built-in GrB_PLUS_FP64 operator, and rdiv is z=y/x in double.
#include "GB_mex.h"
#define USAGE "C = GB_mex_rdiv (A, B, axb_method, cprint)"
#define FREE_ALL \
{ \
GrB_Matrix_free_(&A) ; \
GrB_Matrix_free_(&B) ; \
GrB_Matrix_free_(&C) ; \
GrB_BinaryOp_free_(&My_rdiv) ; \
GrB_Semiring_free_(&My_plus_rdiv) ; \
GB_mx_put_global (true) ; \
}
//------------------------------------------------------------------------------
GrB_Info info ;
bool malloc_debug = false ;
bool ignore = false, ignore1 = false, ignore2 = false ;
bool cprint = false ;
GrB_Matrix A = NULL, B = NULL, C = NULL, MT = NULL ;
int64_t anrows = 0 ;
int64_t ancols = 0 ;
int64_t bnrows = 0 ;
int64_t bncols = 0 ;
GrB_Desc_Value AxB_method = GxB_DEFAULT ;
struct GB_Matrix_opaque MT_header, C_header ;
GrB_Info axb (GB_Context Context, bool cprint) ;
GrB_Semiring My_plus_rdiv = NULL ;
GrB_BinaryOp My_rdiv = NULL ;
void my_rdiv (double *z, const double *x, const double *y) ;
void my_rdiv (double *z, const double *x, const double *y)
{
(*z) = (*y) / (*x) ;
}
#define MY_RDIV \
"void my_rdiv (double *z, const double *x, const double *y)\n" \
"{\n" \
" (*z) = (*y) / (*x) ;\n" \
"}"
//------------------------------------------------------------------------------
GrB_Info axb (GB_Context Context, bool cprint)
{
// create the rdiv operator
// info = GrB_BinaryOp_new (&My_rdiv,
// (GxB_binary_function) my_rdiv, GrB_FP64, GrB_FP64, GrB_FP64) ;
info = GxB_BinaryOp_new (&My_rdiv,
(GxB_binary_function) my_rdiv, GrB_FP64, GrB_FP64, GrB_FP64,
"my_rdiv", MY_RDIV) ;
if (info != GrB_SUCCESS) return (info) ;
GrB_BinaryOp_wait_(My_rdiv, GrB_MATERIALIZE) ;
if (info != GrB_SUCCESS) return (info) ;
info = GrB_Semiring_new (&My_plus_rdiv, GxB_PLUS_FP64_MONOID, My_rdiv) ;
if (info != GrB_SUCCESS)
{
GrB_BinaryOp_free_(&My_rdiv) ;
return (info) ;
}
MT = GB_clear_static_header (&MT_header) ;
C = GB_clear_static_header (&C_header) ;
// C = A*B
info = GB_AxB_meta (C, NULL, // C cannot be computed in place
false, // C_replace
true, // CSC
MT, // no MT returned
&ignore1, // M_transposed will be false
NULL, // no Mask
false, // mask not complemented
false, // mask not structural
NULL, // no accum
A, B,
My_plus_rdiv,
false, // no A transpose
false, // no B transpose
false, // no flipxy
&ignore, // mask_applied
&ignore2, // done_in_place
AxB_method,
true, // do the sort
Context) ;
if (C != NULL)
{
if (cprint) GxB_Matrix_fprint_(C, GxB_COMPLETE, NULL) ;
}
GrB_BinaryOp_free_(&My_rdiv) ;
GrB_Semiring_free_(&My_plus_rdiv) ;
return (info) ;
}
//------------------------------------------------------------------------------
void mexFunction
(
int nargout,
mxArray *pargout [ ],
int nargin,
const mxArray *pargin [ ]
)
{
info = GrB_SUCCESS ;
malloc_debug = GB_mx_get_global (true) ;
ignore = false ;
ignore1 = false ;
ignore2 = false ;
A = NULL ;
B = NULL ;
C = NULL ;
My_rdiv = NULL ;
My_plus_rdiv = NULL ;
GB_CONTEXT (USAGE) ;
// check inputs
if (nargout > 1 || nargin < 2 || nargin > 4)
{
mexErrMsgTxt ("Usage: " USAGE) ;
}
#define GET_DEEP_COPY ;
#define FREE_DEEP_COPY ;
// get A and B
A = GB_mx_mxArray_to_Matrix (pargin [0], "A", false, true) ;
B = GB_mx_mxArray_to_Matrix (pargin [1], "B", false, true) ;
if (A == NULL || B == NULL)
{
FREE_ALL ;
mexErrMsgTxt ("failed") ;
}
if (!A->is_csc || !B->is_csc)
{
mexErrMsgTxt ("A and B must be in CSC format") ;
}
// get the axb_method
// 0 or not present: default
// 1001: Gustavson
// 1003: dot
// 1004: hash
// 1005: saxpy
GET_SCALAR (2, GrB_Desc_Value, AxB_method, GxB_DEFAULT) ;
// get the cprint flag
GET_SCALAR (3, bool, cprint, false) ;
if (! ((AxB_method == GxB_DEFAULT) ||
(AxB_method == GxB_AxB_GUSTAVSON) ||
(AxB_method == GxB_AxB_HASH) ||
(AxB_method == GxB_AxB_SAXPY) ||
(AxB_method == GxB_AxB_DOT)))
{
mexErrMsgTxt ("unknown method") ;
}
// determine the dimensions
anrows = GB_NROWS (A) ;
ancols = GB_NCOLS (A) ;
bnrows = GB_NROWS (B) ;
bncols = GB_NCOLS (B) ;
if (ancols != bnrows)
{
FREE_ALL ;
mexErrMsgTxt ("invalid dimensions") ;
}
METHOD (axb (Context, cprint)) ;
// return C
pargout [0] = GB_mx_Matrix_to_mxArray (&C, "C AxB result", false) ;
FREE_ALL ;
}
|
