File: GB_iso_AxB.c

package info (click to toggle)
suitesparse-graphblas 7.4.0%2Bdfsg-2
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid, trixie
  • size: 67,112 kB
  • sloc: ansic: 1,072,243; cpp: 8,081; sh: 512; makefile: 506; asm: 369; python: 125; awk: 10
file content (324 lines) | stat: -rw-r--r-- 12,338 bytes parent folder | download | duplicates (3)
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
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
//------------------------------------------------------------------------------
// GB_iso_AxB: check for iso result for C=A*B and compute the iso scalar for C
//------------------------------------------------------------------------------

// SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2022, All Rights Reserved.
// SPDX-License-Identifier: Apache-2.0

//------------------------------------------------------------------------------

// Return true if C=A*B results in an iso matrix C, and return the iso value of
// C.  The type of the matrix C and scalar c is semiring->add->ztype.

// If both A and B are full and iso, then C is also full and iso, for nearly
// all semirings.  The inner dimension of the matrix multiply is required to
// compute the iso value of C.  Assuming all matrices are CSC:

//  C = A*B     n = A->vdim == B->vlen
//  C = A'*B    n = A->vlen == B->vlen
//  C = A*B'    n = A->vdim == B->vdim
//  C = A'*B'   n = A->vlen == B->vdim

#include "GB_mxm.h"
#include "GB_reduce.h"
#include "GB_binop.h"

//------------------------------------------------------------------------------
// GB_iso_mult: c = mult(a,b) or c = mult(b,a)
//------------------------------------------------------------------------------

static void GB_iso_mult         // c = mult(a,b) or c=mult(b,a)
(
    GB_void *restrict c,        // c has type zcode, and size zsize
    const GB_void *restrict a, const GB_Type_code acode, const size_t asize,
    const GB_void *restrict b, const GB_Type_code bcode, const size_t bsize,
    GxB_binary_function fmult,
    const bool flipxy,
    const GB_Type_code xcode, const size_t xsize,
    const GB_Type_code ycode, const size_t ysize,
    const GB_Type_code zcode, const size_t zsize
)
{
    if (flipxy)
    { 
        // c = mult(b,a)
        GB_iso_mult (c, b, bcode, bsize, a, acode, asize, fmult, false,
            xcode, xsize, ycode, ysize, zcode, zsize) ;
    }
    else
    {
        if (fmult == NULL)
        { 
            // fmult is the implicit FIRST operator from GB_reduce_to_vector
            // c = (ztype) a
            GB_cast_scalar (c, zcode, a, acode, asize) ;
        }
        else if (acode == xcode && bcode == ycode)
        { 
            // c = fmult (a,b)
            fmult (c, a, b) ;
        }
        else
        { 
            // x = (xtype) a
            GB_void x [GB_VLA(xsize)] ;
            GB_cast_scalar (x, xcode, a, acode, asize) ;
            // y = (ytype) b
            GB_void y [GB_VLA(ysize)] ;
            GB_cast_scalar (y, ycode, b, bcode, bsize) ;
            // c = fmult (x,y)
            fmult (c, x, y) ;
        }
    }
}

//------------------------------------------------------------------------------
// GB_iso_AxB
//------------------------------------------------------------------------------

bool GB_iso_AxB             // C = A*B, return true if C is iso
(
    // output
    GB_void *restrict c,    // output scalar of iso array (not computed if NULL)
    // input
    GrB_Matrix A,           // input matrix
    GrB_Matrix B,           // input matrix
    uint64_t n,             // inner dimension of the matrix multiply
    GrB_Semiring semiring,  // semiring
    bool flipxy,            // true if z=fmult(b,a), false if z=fmult(a,b)
    bool ignore_monoid      // rowscale and colscale do not use the monoid
)
{

    //--------------------------------------------------------------------------
    // get inputs
    //--------------------------------------------------------------------------

    ASSERT_MATRIX_OK (A, "A for GB_iso_AxB", GB0) ;
    ASSERT_MATRIX_OK (B, "B for GB_iso_AxB", GB0) ;
    ASSERT_SEMIRING_OK (semiring, "semiring for GB_iso_AxB", GB0) ;

    //--------------------------------------------------------------------------
    // quick return if multop is positional
    //--------------------------------------------------------------------------

    GB_Opcode add_binop_code = semiring->add->op->opcode ;
    const GrB_BinaryOp multiply = semiring->multiply ;

    if (GB_OP_IS_POSITIONAL (multiply))
    { 
        // C is not iso if the multiply op is positional
        return (false) ;
    }

    //--------------------------------------------------------------------------
    // get the binary operator and the types of C, A, and B
    //--------------------------------------------------------------------------

    const GxB_binary_function fmult = multiply->binop_function ;
    GB_Opcode mult_binop_code = multiply->opcode ;
    ASSERT (GB_IMPLIES (fmult == NULL, mult_binop_code == GB_FIRST_binop_code));

    const GrB_Type xtype = multiply->xtype ;
    const GrB_Type ytype = multiply->ytype ;
    const GrB_Type ztype = multiply->ztype ;

    const GB_Type_code xcode = xtype->code ;
    const GB_Type_code ycode = ytype->code ;
    const GB_Type_code zcode = ztype->code ;
    const GB_Type_code acode = A->type->code ;
    const GB_Type_code bcode = B->type->code ;

    const size_t xsize = xtype->size ;
    const size_t ysize = ytype->size ;
    const size_t zsize = ztype->size ;
    const size_t asize = A->type->size ;
    const size_t bsize = B->type->size ;

    if (zcode == GB_BOOL_code)
    { 
        // rename a boolean monoid:
        // MIN_BOOL and TIMES_BOOL monoids become LAND
        // MAX_BOOL and PLUS_BOOL monoids become LOR
        add_binop_code = GB_boolean_rename (add_binop_code) ;
    }

    if (xcode == GB_BOOL_code)
    { 
        // rename a boolean multiply op:
        // DIV becomes FIRST, RDIV becomes SECOND; all other renaming has no
        // effect on this method. 
        mult_binop_code = GB_boolean_rename (mult_binop_code) ;
    }

    // "nice" monoids have the property that reducing a set of iso values to a
    // single result doesn't change the result: ANY, LAND, LOR, BAND, BOR, MIN
    // and MAX.  That is, x == reduce ([x x x x x ... x x x x]) holds for all
    // these monoids.  The monoids that do not fall into this "nice" category
    // are PLUS, TIMES, EQ (LXNOR), LXOR, BXOR, and BXNOR.  For row/col scaling,
    // all monoids are "nice" since they aren't used.
    const bool nice_monoid = ignore_monoid ||
        add_binop_code == GB_ANY_binop_code  ||
        add_binop_code == GB_LAND_binop_code ||
        add_binop_code == GB_LOR_binop_code  ||
        add_binop_code == GB_BAND_binop_code ||
        add_binop_code == GB_BOR_binop_code  ||
        add_binop_code == GB_MAX_binop_code  ||
        add_binop_code == GB_MIN_binop_code ;

    // Nearly all cases where C is iso require a "nice" monoid, with the
    // exception of the EQ_PAIR and TIMES_PAIR semirings, which are the same
    // as ANY_PAIR semirings.
    const bool nice_with_pair = nice_monoid ||
        add_binop_code == GB_EQ_binop_code  ||
        add_binop_code == GB_TIMES_binop_code ;

    // the FIRST or ANY multiply ops can both produce a FIRST result
    const bool first = (mult_binop_code == GB_ANY_binop_code) ||
        (mult_binop_code ==
            (flipxy ? GB_SECOND_binop_code : GB_FIRST_binop_code)) ;

    // the SECOND or ANY multiply ops can both produce a SECOND result
    const bool second = (mult_binop_code == GB_ANY_binop_code) ||
        (mult_binop_code ==
            (flipxy ? GB_FIRST_binop_code : GB_SECOND_binop_code)) ;

    //--------------------------------------------------------------------------
    // determine if C is iso
    //--------------------------------------------------------------------------

    // A and B are treated as if iso if they have 1 entry and are not bitmap
    const bool A_iso = A->iso || (GB_nnz (A) == 1 && !GB_IS_BITMAP (A)) ;
    const bool B_iso = B->iso || (GB_nnz (B) == 1 && !GB_IS_BITMAP (B)) ;

    if (nice_with_pair && mult_binop_code == GB_PAIR_binop_code)
    {

        //----------------------------------------------------------------------
        // C is iso, with c = 1
        //----------------------------------------------------------------------

        if (c != NULL)
        { 
            GB_cast_one (c, zcode) ;
        }
        return (true) ;

    }
    else if (B_iso && nice_monoid && second)
    {

        //----------------------------------------------------------------------
        // C is iso, with c = b
        //----------------------------------------------------------------------

        if (c != NULL)
        { 
            if (zcode == ycode && bcode == ycode)
            { 
                // c = Bx [0]
                memcpy (c, B->x, zsize) ;
            }
            else
            { 
                // c = (ztype) ((ytype) Bx [0])
                GB_void y [GB_VLA(ysize)] ;
                GB_cast_scalar (y, ycode, B->x, bcode, bsize) ;
                GB_cast_scalar (c, zcode, y, ycode, ysize) ;
            }
        }
        return (true) ;

    }
    else if (A_iso && nice_monoid && first)
    {

        //----------------------------------------------------------------------
        // C is iso, with c = a
        //----------------------------------------------------------------------

        if (c != NULL)
        { 
            if (zcode == xcode && acode == xcode)
            { 
                // c = Ax [0]
                memcpy (c, A->x, zsize) ;
            }
            else
            { 
                // c = (ztype) ((xtype) Ax [0])
                GB_void x [GB_VLA(xsize)] ;
                GB_cast_scalar (x, xcode, A->x, acode, asize) ;
                GB_cast_scalar (c, zcode, x, xcode, xsize) ;
            }
        }
        return (true) ;

    }
    else if (A_iso && B_iso)
    {

        //----------------------------------------------------------------------
        // both A and B are iso
        //----------------------------------------------------------------------

        GB_void *Ax = (GB_void *) A->x ;
        GB_void *Bx = (GB_void *) B->x ;
        if (nice_monoid)
        {

            //------------------------------------------------------------------
            // C is iso, with c = fmult(a,b), for any fmult, incl. user-defined
            //------------------------------------------------------------------

            if (c != NULL)
            { 
                GB_iso_mult (c, Ax, acode, asize, Bx, bcode, bsize,
                    fmult, flipxy, xcode, xsize, ycode, ysize, zcode, zsize) ;
            }
            return (true) ;

        }
        else if (GB_as_if_full (A) && GB_as_if_full (B))
        {

            //------------------------------------------------------------------
            // C = A*B where A and B are both full and iso
            //------------------------------------------------------------------

            // If A and B are both full and iso, then C is also full and iso,
            // for any semiring (including user-defined) except those with a
            // positional multiplicative operator.  Each entry C(i,j) is the
            // reduction of n copies of the single iso scalar t, where t =
            // A(i,k)*B(k,j) is iso-valued for any i, j, or k, assuming n is
            // the inner dimension of the C=A*B matrix multiply.

            if (c != NULL)
            { 
                // t = A(i,k)*B(k,j)
                GB_void t [GB_VLA(zsize)] ;
                GB_iso_mult (t, Ax, acode, asize, Bx, bcode, bsize,
                    fmult, flipxy, xcode, xsize, ycode, ysize, zcode, zsize) ;

                // reduce n copies of t to the single scalar c, in O(log(n))
                GxB_binary_function freduce = semiring->add->op->binop_function;
                GB_iso_reduce_worker (c, freduce, t, n, zsize) ;
            }

            // the total time to compute C=A*B where all matrices are n-by-n
            // and full is thus O(log(n)), much smaller than O(n^3) for the
            // conventional matrix-multiply algorithm.  It would be possible to
            // reduce the time still further, since most reductions of n copies
            // of t can be done in O(1) time, but the O(log(n)) method works
            // for any monoid, including user-defined ones.
            return (true) ;
        }
    }

    //--------------------------------------------------------------------------
    // otherwise, C is not iso
    //--------------------------------------------------------------------------

    return (false) ;
}