1function C = GB_spec_op (op, A, B)
2%GB_SPEC_OP apply a unary or binary operator
3%
4% Apply a binary operator z = f (x,y) element-wise to x and y, or a unary
5% operator z = f(x) just x.  The operator op is any built-in GraphBLAS
6% operator.
7%
8% op or op.opname is a string with just the operator name.  Valid names of
9% binary operators are 'first', 'second', 'min', 'max', 'plus', 'minus',
10% 'rminus', 'times', 'div', 'rdiv', 'eq', 'ne', 'gt', 'lt', 'ge', 'le', 'or',
11% 'and', 'xor'.  'iseq', 'isne', 'isgt', 'islt', 'isge', 'le', 'pair', 'any',
12% 'pow', ('bitget' or 'bget'), ('bitset' or 'bset'), ('bitclr' or 'bclr'),
13% ('bitand' or 'band'), ('bitor' or 'bor'), ('bitxor' or 'bxor'), ('bitxnor',
14% 'bxnor'), ('bitshift' or 'bshift'), ('bitnot' or 'bitcmp'), 'atan2', 'hypot',
15% ('ldexp' or 'pow2'), ('complex', 'cmplx').
16%
17% Unary operators are 'one', 'identity', 'ainv', 'abs', 'minv', 'not', 'bnot',
18% 'sqrt', 'log', 'exp', 'sin', 'cos', 'tan', 'asin', 'acos', 'atan', 'sinh',
19% 'cosh', 'tanh', 'asinh', 'acosh', 'atanh', 'signum', 'ceil', 'floor',
20% 'round', ('trunc' or 'fix'), 'exp2', 'expm1', 'log10', 'log2', ('lgamma' or
21% 'gammaln'), ('tgamma' or 'gamma'), 'erf', 'erfc', 'frexpx', 'frexpe', 'conj',
22% ('creal' or 'real'), ('cimag' or 'imag'), ('carg' or 'angle'), 'isinf',
23% 'isnan', 'isfinite'.
24%
25% op.optype: 'logical', 'int8', 'uint8', 'int16', 'uint16', 'int32', 'uint32',
26% 'int64', 'uint64', 'single', 'double', 'single complex' or 'double complex'.
27%
28% The class of z is the same as the class of the output of the operator, which
29% is op.optype except for: (1) 'eq', 'ne', 'gt', 'lt', 'ge', 'le', in which
30% case z is logical, (2) 'complex', where x and y are real and z is complex,
31% (3) bitshift (where x has the optype and y is int8).
32%
33% Intrinsic MATLAB operators are used as much as possible, so as to test
34% GraphBLAS operators.  Some must be done in GraphBLAS because the
35% divide-by-zero and overflow rules for integers differs between MATLAB and C.
36% Also, typecasting in MATLAB and GraphBLAS differs with underflow and overflow
37% conditions.
38%
39% Positional ops are not computed by this function.
40
41% SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2021, All Rights Reserved.
42% SPDX-License-Identifier: Apache-2.0
43
44% get the operator name and class
45[opname optype ztype xtype ytype] = GB_spec_operator (op, GB_spec_type (A)) ;
46
47if (GB_spec_is_positional (opname))
48    error ('positional op not supported by this funciton') ;
49end
50
51% cast the inputs A and B to the inputs of the operator
52if (~isequal (GB_spec_type (A), xtype))
53    x = GB_mex_cast (A, xtype) ;
54else
55    x = A ;
56end
57
58use_matlab = (isa (x, 'float') && ...
59    (contains (optype, 'single') || contains (optype, 'double'))) ;
60
61if (nargin > 2 && ~ischar (B))
62    if (~isequal (GB_spec_type (B), ytype))
63        y = GB_mex_cast (B, ytype) ;
64    else
65        y = B ;
66    end
67    use_matlab = use_matlab && isa (y, 'float') ;
68end
69
70switch opname
71
72    % binary operators, result is ztype
73    case 'first'
74        z = x ;
75    case 'second'
76        z = y ;
77    case 'any'
78        z = y ;
79    case 'pair'
80        z = GB_spec_ones (size (x), ztype) ;
81    case 'min'
82        % min(x,y) in SuiteSparse:GraphBLAS is min(x,y,'omitnan') in MATLAB.
83        % see discussion in SuiteSparse/GraphBLAS/Source/GB.h
84        % z = min (x,y,'omitnan') ;
85        z = GB_mex_op (op, x, y) ;
86    case 'max'
87        % z = max (x,y,'omitnan') ;
88        z = GB_mex_op (op, x, y) ;
89    case 'plus'
90        if (use_matlab)
91            z = x + y ;
92        else
93            z = GB_mex_op (op, x, y) ;
94        end
95    case 'minus'
96        if (use_matlab)
97            z = x - y ;
98        else
99            z = GB_mex_op (op, x, y) ;
100        end
101    case 'rminus'
102        if (use_matlab)
103            z = y - x ;
104        else
105            z = GB_mex_op (op, x, y) ;
106        end
107    case 'times'
108        if (use_matlab)
109            z = x .* y ;
110        else
111            z = GB_mex_op (op, x, y) ;
112        end
113    case 'div'
114        if (use_matlab)
115            z = x ./ y ;
116        else
117            z = GB_mex_op (op, x, y) ;
118        end
119    case 'rdiv'
120        if (use_matlab)
121            z = y ./ x ;
122        else
123            z = GB_mex_op (op, x, y) ;
124        end
125    case 'pow'
126        if (use_matlab)
127            z = x .^ y ;
128        else
129            z = GB_mex_op (op, x, y) ;
130        end
131
132    % 6 binary comparison operators (result is ztype)
133    case 'iseq'
134        z = GB_mex_cast (x == y, ztype) ;
135    case 'isne'
136        z = GB_mex_cast (x ~= y, ztype) ;
137    case 'isgt'
138        z = GB_mex_cast (x >  y, ztype) ;
139    case 'islt'
140        z = GB_mex_cast (x <  y, ztype) ;
141    case 'isge'
142        z = GB_mex_cast (x >= y, ztype) ;
143    case 'isle'
144        z = GB_mex_cast (x <= y, ztype) ;
145
146    % 6 binary comparison operators (result is boolean)
147    case 'eq'
148        z = (x == y) ;
149    case 'ne'
150        z = (x ~= y) ;
151    case 'gt'
152        z = (x >  y) ;
153    case 'lt'
154        z = (x <  y) ;
155    case 'ge'
156        z = (x >= y) ;
157    case 'le'
158        z = (x <= y) ;
159
160    % 3 binary logical operators (result is ztype)
161    case 'or'
162        z = GB_mex_cast ((x ~= 0) | (y ~= 0), ztype) ;
163    case 'and'
164        z = GB_mex_cast ((x ~= 0) & (y ~= 0), ztype) ;
165    case 'xor'
166        z = GB_mex_cast ((x ~= 0) ~= (y ~= 0), ztype) ;
167
168    % bitwise operators
169    case { 'bitget', 'bget' }
170        bits = GB_spec_nbits (ztype) ;
171        m = (y > 0 & y <= bits) ;
172        t = bitget (x (m), y (m), ztype) ;
173        z = zeros (size (x), ztype) ;
174        z (m) = t ;
175
176    case { 'bitset', 'bset' }
177        bits = GB_spec_nbits (ztype) ;
178        m = (y > 0 & y <= bits) ;
179        t = bitset (x (m), y (m), 1, ztype) ;
180        z = x ;
181        z (m) = t ;
182
183    case { 'bitclr', 'bclr' }
184        bits = GB_spec_nbits (ztype) ;
185        m = (y > 0 & y <= bits) ;
186        t = bitset (x (m), y (m), 0, ztype) ;
187        z = x ;
188        z (m) = t ;
189
190    case { 'bitand', 'band' }
191        z = bitand (x, y, ztype) ;
192    case { 'bitor', 'bor' }
193        z = bitor (x, y, ztype) ;
194    case { 'bitxor', 'bxor' }
195        z = bitxor (x, y, ztype) ;
196    case { 'bitxnor', 'bxnor' }
197        z = bitcmp (bitxor (x, y, ztype), ztype) ;
198    case { 'bitshift', 'bshift' }
199        z = bitshift (x, y, ztype) ;
200    case { 'bitnot', 'bitcmp', 'bnot', 'bcmp' }
201        z = bitcmp (x, ztype) ;
202
203    case 'atan2'
204        z = atan2 (x,y) ;
205
206    case 'hypot'
207        z = hypot (x,y) ;
208
209    case { 'ldexp', 'pow2' }
210        z = pow2 (x,y) ;
211
212    case { 'fmod', 'rem' }
213        % see ANSI C11 fmod function
214        % the MATLAB rem differs slightly from the ANSI C11 fmod,
215        % if x/y is O(eps) smaller than an integer.
216        z = rem (x,y) ;
217
218    case { 'remainder' }
219        % see ANSI C11 remainder function
220        m = (y ~= 0 & x ~= y)  ;
221        z = nan (size (x), ztype) ;
222        z (x == y) = 0 ;
223        z (m) = x (m) - round (x (m) ./ y (m)) .* y (m) ;
224
225    case { 'copysign' }
226        % see ANSI C11 copysign function
227        z = abs (x) .* (2 * double (y >= 0) - 1) ;
228
229    case { 'complex', 'cmplx' }
230        z = complex (x,y) ;
231
232    % unary operators (result is ztype)
233    case 'one'
234        z = GB_mex_cast (1, ztype) ;
235    case 'identity'
236        z = x ;
237    case 'ainv'
238        if (use_matlab)
239            z = -x ;
240        else
241            z = GB_mex_op (op, x) ;
242        end
243    case 'abs'
244        if (use_matlab)
245            z = abs (x) ;
246        else
247            z = GB_mex_op (op, x) ;
248        end
249    case 'minv'
250        if (use_matlab)
251            z = 1 ./ x ;
252        else
253            z = GB_mex_op (op, x) ;
254        end
255    case 'not'
256        z = GB_mex_cast (~(x ~= 0), ztype) ;
257
258    case 'bnot'
259        z = bitcmp (x) ;
260
261    case 'sqrt'
262        z = sqrt (x) ;
263
264    case 'log'
265        z = log (x) ;
266
267    case 'exp'
268        z = exp (x) ;
269
270    case 'sin'
271        z = sin (x) ;
272
273    case 'cos'
274        z = cos (x) ;
275
276    case 'tan'
277        z = tan (x) ;
278
279    case 'asin'
280        z = asin (x) ;
281
282    case 'acos'
283        z = acos (x) ;
284
285    case 'atan'
286        z = atan (x) ;
287
288    case 'sinh'
289        z = sinh (x) ;
290
291    case 'cosh'
292        z = cosh (x) ;
293
294    case 'tanh'
295        z = tanh (x) ;
296
297    case 'asinh'
298        z = asinh (x) ;
299
300    case 'acosh'
301        z = acosh (x) ;
302
303    case 'atanh'
304        z = atanh (x) ;
305
306    case { 'sign', 'signum' }
307        z = sign (x) ;
308
309    case 'ceil'
310        z = ceil (x) ;
311
312    case 'floor'
313        z = floor (x) ;
314
315    case 'round'
316        z = round (x) ;
317
318    case { 'trunc', 'fix' }
319        z = fix (x) ;
320
321    case { 'exp2' }
322        z = 2.^x ;
323
324    case 'expm1'
325        z = expm1 (x) ;
326
327    case 'log10'
328        z = log10 (x) ;
329
330    case 'log2'
331        z = log2 (x) ;
332
333    case 'log1p'
334        z = log1p (x) ;
335
336    case { 'lgamma', 'gammaln' }
337        z = gammaln (x) ;
338
339    case { 'tgamma', 'gamma' }
340        z = gamma (x) ;
341
342    case 'erf'
343        z = erf (x) ;
344
345    case 'erfc'
346        z = erfc (x) ;
347
348    case 'frexpx'
349        [z,~] = log2 (x) ;
350
351    case 'frexpe'
352        [~,z] = log2 (x) ;
353
354    case 'conj'
355        z = conj (x) ;
356
357    case { 'creal', 'real' }
358        z = real (x) ;
359
360    case { 'cimag', 'imag' }
361        z = imag (x) ;
362
363    case { 'carg', 'angle' }
364        z = angle (x) ;
365
366    case 'isinf'
367        z = isinf (x) ;
368
369    case 'isnan'
370        z = isnan (x) ;
371
372    case 'isfinite'
373        z = isfinite (x) ;
374
375    otherwise
376        opname
377        error ('unknown op') ;
378end
379
380if (~isequal (ztype, GB_spec_type (z)))
381    z = GB_mex_cast (z, ztype) ;
382end
383
384C = z ;
385
386