Skip to content

Commit

Permalink
using GrB_UINT32 in FastSV if n <= INT32_MAX
Browse files Browse the repository at this point in the history
  • Loading branch information
@DrTimothyAldenDavis
DrTimothyAldenDavis committed Nov 27, 2021
1 parent 03cae93 commit 0b0cb0a
Show file tree
Hide file tree
Showing 2 changed files with 131 additions and 89 deletions.
206 changes: 117 additions & 89 deletions src/algorithm/LG_CC_FastSV6.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,74 +37,77 @@
#endif

//------------------------------------------------------------------------------
// Reduce_assign: w (Px) += s, using MIN as the "+=" accum operator
// Reduce_assign: parent (Px) += mngp, using MIN as the "+=" accum operator
//------------------------------------------------------------------------------

// The Px array of size n is the non-opaque version of the parent vector, where
// The Px array of size n is the non-opaque copy of the parent vector, where
// i = Px [j] if the parent of node j is node i. It can thus have duplicates.
// The vectors w and s are full (all entries present). This function computes
// the following, which is done explicitly in the Reduce_assign function in
// LG_CC_Boruvka:
// The vectors parent and mngp are full (all entries present). This function
// computes the following, which is done explicitly in the Reduce_assign
// function in LG_CC_Boruvka:
//
// for (j = 0 ; j < n ; j++)
// {
// uint64_t i = Px [j] ;
// w [i] = min (w [i], s [j]) ;
// parent [i] = min (parent [i], mngp [j]) ;
// }
//
// If C(i,j) = true where i == Px [j], then this can be written as:
// If C(i,j) is present where i == Px [j], then this can be written as:
//
// w = min (w, C*s)
// parent = min (parent, C*mngp)
//
// when using the min_second semiring. This can be done efficiently where
// when using the min_2nd semiring. This can be done efficiently where
// because C can be constructed in O(1) time and O(1) additional space (not
// counting the prior Cp, Ci, and Cx arrays), when using the SuiteSparse
// pack/unpack move constructors.
// pack/unpack move constructors. The min_2nd semiring ignores the values of
// C and operates only on the structure, so its values are not relevant.
// C is thus chosen as a GrB_BOOL array where Cx [0] = false.

static inline GrB_Info Reduce_assign
(
GrB_Vector w, // vector of size n, all entries present
GrB_Vector s, // vector of size n, all entries present
// input/output:
GrB_Vector parent, // vector of size n, all entries present
// input:
GrB_BinaryOp min, // min operator (uint32 or uint64)
GrB_Semiring min_2nd, // min_2nd semiring (uint32 or uint64)
GrB_Vector mngp, // vector of size n, all entries present
GrB_Matrix C, // boolean matrix of size n-by-n
GrB_Index **Cp_handle, // array of size n+1, equal to 0:n
GrB_Index **Ci_handle, // Px array of size n, can have duplicates
bool **Cx_handle, // array of size 1, equal to true
GrB_Index **Ci_handle, // Px array of size n, always uint64
void **Cx_handle, // array of size 1, contents not accessed
char *msg
)
{
// size of Cp, Ci, and Cx in bytes
GrB_Index n ;
GrB_TRY (GrB_Vector_size (&n, w)) ;
GrB_TRY (GrB_Vector_size (&n, parent)) ;
GrB_Index Cp_size = (n+1) * sizeof (GrB_Index) ;
GrB_Index Ci_size = n * sizeof (GrB_Index) ;
GrB_Index Cx_size = sizeof (bool) ;

// pack Cp, Ci, and Cx into a matrix C with C(i,j) = true if Ci(j) == i
bool iso = true ;
bool jumbled = false ;
GrB_TRY (GxB_Matrix_pack_CSC (C, Cp_handle, Ci_handle, (void **) Cx_handle,
// pack Cp, Ci, and Cx into a matrix C with C(i,j) present if Ci(j) == i
bool iso = true, jumbled = false ;
GrB_TRY (GxB_Matrix_pack_CSC (C, Cp_handle, Ci_handle, Cx_handle,
Cp_size, Ci_size, Cx_size, iso, jumbled, NULL)) ;

// w = min (w, C*s) using the MIN_SECOND semiring
GrB_TRY (GrB_mxv (w, NULL, GrB_MIN_UINT64,
GrB_MIN_SECOND_SEMIRING_UINT64, C, s, NULL)) ;
// parent = min (parent, C*mngp) using the MIN_SECOND semiring
GrB_TRY (GrB_mxv (parent, NULL, min, min_2nd, C, mngp, NULL)) ;

// unpack the contents of C
GrB_TRY (GxB_Matrix_unpack_CSC (C, Cp_handle, Ci_handle, (void **)Cx_handle,
GrB_TRY (GxB_Matrix_unpack_CSC (C, Cp_handle, Ci_handle, Cx_handle,
&Cp_size, &Ci_size, &Cx_size, &iso, &jumbled, NULL)) ;

return (GrB_SUCCESS) ; // yay! It works!
return (GrB_SUCCESS) ;
}

//------------------------------------------------------------------------------
// LG_CC_FastSV6
//------------------------------------------------------------------------------

// The output of LG_CC_FastSV* is a vector component, where
// component(i)=s if node i is in the connected compononent whose
// representative node is node s. If s is a representative, then
// component(s)=s. The number of connected components in the graph G is the
// number of representatives.
// The output of LG_CC_FastSV* is a vector component, where component(i)=s if
// node i is in the connected compononent whose representative is node s. If s
// is a representative, then component(s)=s. The number of connected
// components in the graph G is the number of representatives.

#undef LAGraph_FREE_WORK
#define LAGraph_FREE_WORK \
Expand All @@ -113,8 +116,8 @@ static inline GrB_Info Reduce_assign
LAGraph_Free ((void **) &Tj) ; \
LAGraph_Free ((void **) &Tx) ; \
LAGraph_Free ((void **) &Cp) ; \
LAGraph_Free ((void **) &Cx) ; \
LAGraph_Free ((void **) &Px) ; \
LAGraph_Free ((void **) &Cx) ; \
LAGraph_Free ((void **) &ht_key) ; \
LAGraph_Free ((void **) &ht_count) ; \
LAGraph_Free ((void **) &count) ; \
Expand All @@ -125,7 +128,6 @@ static inline GrB_Info Reduce_assign
GrB_free (&gp) ; \
GrB_free (&mngp) ; \
GrB_free (&gp_new) ; \
GrB_free (&mod) ; \
}

#undef LAGraph_FREE_ALL
Expand Down Expand Up @@ -157,14 +159,14 @@ int LG_CC_FastSV6 // SuiteSparse:GraphBLAS method, with GxB extensions

LG_CLEAR_MSG ;

int64_t *ht_key = NULL, *ht_count = NULL, *range = NULL ;
GrB_Index n, nnz, Cp_size = 0,
*Px = NULL, *Cp = NULL, *count = NULL, *Tp = NULL, *Tj = NULL ;
GrB_Vector parent = NULL, gp_new = NULL, mngp = NULL, mod = NULL, gp = NULL,
t = NULL, y = NULL ;
int64_t *range = NULL ;
GrB_Index n, nnz, Cp_size = 0, *ht_key = NULL, *Px = NULL, *Cp = NULL,
*count = NULL, *Tp = NULL, *Tj = NULL ;
GrB_Vector parent = NULL, gp_new = NULL, mngp = NULL, gp = NULL, t = NULL,
y = NULL ;
GrB_Matrix T = NULL, C = NULL ;
bool *Cx = NULL ;
void *Tx = NULL ;
void *Tx = NULL, *Cx = NULL ;
int *ht_count = NULL ;

LG_CHECK (LAGraph_CheckGraph (G, msg), GrB_INVALID_OBJECT,
"graph is invalid") ;
Expand All @@ -191,6 +193,36 @@ int LG_CC_FastSV6 // SuiteSparse:GraphBLAS method, with GxB extensions
GrB_TRY (GrB_Matrix_nrows (&n, A)) ;
GrB_TRY (GrB_Matrix_nvals (&nnz, A)) ;

// determine the integer type, operators, and semirings to use
GrB_Type Uint, Int ;
GrB_IndexUnaryOp ramp ;
GrB_Semiring min_2nd, min_2ndi ;
GrB_BinaryOp min, ne, imin ;
if (n > INT32_MAX)
{
// use 64-bit integers throughout
Uint = GrB_UINT64 ;
Int = GrB_INT64 ;
ramp = GrB_ROWINDEX_INT64 ;
min = GrB_MIN_UINT64 ;
imin = GrB_MIN_INT64 ;
ne = GrB_NE_UINT64 ;
min_2nd = GrB_MIN_SECOND_SEMIRING_UINT64 ;
min_2ndi = GxB_MIN_SECONDI_INT64 ;
}
else
{
// use 32-bit integers, except for Px and for constructing the matrix C
Uint = GrB_UINT32 ;
Int = GrB_INT32 ;
ramp = GrB_ROWINDEX_INT32 ;
min = GrB_MIN_UINT32 ;
imin = GrB_MIN_INT32 ;
ne = GrB_NE_UINT32 ;
min_2nd = GrB_MIN_SECOND_SEMIRING_UINT32 ;
min_2ndi = GxB_MIN_SECONDI_INT32 ;
}

// FASTSV_SAMPLES: number of samples to take from each row A(i,:)
#define FASTSV_SAMPLES 4
bool sampling = (n * FASTSV_SAMPLES * 2 < nnz && n > 1024) ;
Expand All @@ -201,15 +233,13 @@ int LG_CC_FastSV6 // SuiteSparse:GraphBLAS method, with GxB extensions
nthreads = LAGraph_MIN (nthreads, n / 16) ;
nthreads = LAGraph_MAX (nthreads, 1) ;

GrB_TRY (GrB_Vector_new (&gp_new, GrB_UINT64, n)) ;
GrB_TRY (GrB_Vector_new (&mod, GrB_BOOL, n)) ;
GrB_TRY (GrB_Vector_new (&gp_new, Uint, n)) ;

Cx = (bool *) LAGraph_Malloc (1, sizeof (bool)) ;
Px = LAGraph_Malloc (n, sizeof (uint64_t)) ;
Cx = (void *) LAGraph_Calloc (1, sizeof (bool)) ;
Px = (GrB_Index *) LAGraph_Malloc (n, sizeof (GrB_Index)) ;
LG_CHECK (Px == NULL || Cx == NULL, GrB_OUT_OF_MEMORY, "out of memory") ;
Cx [0] = true ;

// create Cp = 0:n and the empty C matrix
// create Cp = 0:n (always 64-bit) and the empty C matrix
GrB_TRY (GrB_Matrix_new (&C, GrB_BOOL, n, n)) ;
GrB_TRY (GrB_Vector_new (&t, GrB_INT64, n+1)) ;
GrB_TRY (GrB_assign (t, NULL, NULL, 0, GrB_ALL, n+1, NULL)) ;
Expand All @@ -227,27 +257,27 @@ int LG_CC_FastSV6 // SuiteSparse:GraphBLAS method, with GxB extensions
// is the minimum index j, so only the first entry in A(i,:) needs to be
// considered for each row i.

GrB_TRY (GrB_Vector_new (&t, GrB_INT64, n)) ;
GrB_TRY (GrB_Vector_new (&y, GrB_INT64, n)) ;
GrB_TRY (GrB_Vector_new (&t, Int, n)) ;
GrB_TRY (GrB_Vector_new (&y, Int, n)) ;
GrB_TRY (GrB_assign (t, NULL, NULL, 0, GrB_ALL, n, NULL)) ;
GrB_TRY (GrB_assign (y, NULL, NULL, 0, GrB_ALL, n, NULL)) ;
GrB_TRY (GrB_apply (y, NULL, NULL, GrB_ROWINDEX_INT64, y, 0, NULL)) ;
GrB_TRY (GrB_mxv (y, NULL, GrB_MIN_INT64, GxB_MIN_SECONDI_INT64, A, t,
NULL)) ;
GrB_TRY (GrB_apply (y, NULL, NULL, ramp, y, 0, NULL)) ;
GrB_TRY (GrB_mxv (y, NULL, imin, min_2ndi, A, t, NULL)) ;
GrB_TRY (GrB_free (&t)) ;

// The typecast is required because the ROWINDEX operator and MIN_SECONDI
// do not work in the UINT64 domain, as built-in operators.
// parent = (uint64) y
GrB_TRY (GrB_Vector_new (&parent, GrB_UINT64, n)) ;
// The typecast from Int to Uint is required because the ROWINDEX operator
// and MIN_SECONDI do not work in the UINT* domains, as built-in operators.
// parent = (Uint) y
GrB_TRY (GrB_Vector_new (&parent, Uint, n)) ;
GrB_TRY (GrB_assign (parent, NULL, NULL, y, GrB_ALL, n, NULL)) ;
GrB_TRY (GrB_free (&y)) ;

// copy parent into gp, mngp, and Px. Px is a non-opaque copy of the
// parent GrB_Vector
// copy parent into gp, mngp, and Px. Px is a non-opaque 64-bit copy of the
// parent GrB_Vector. If parent is uint32, GraphBLAS typecasts to uint64.
GrB_TRY (GrB_Vector_extractTuples (NULL, Px, &n, parent)) ;
GrB_TRY (GrB_Vector_dup (&gp, parent)) ;
GrB_TRY (GrB_Vector_dup (&mngp, parent)) ;
GrB_TRY (GrB_Vector_new (&t, GrB_BOOL, n)) ;

//--------------------------------------------------------------------------
// sample phase
Expand All @@ -274,11 +304,11 @@ int LG_CC_FastSV6 // SuiteSparse:GraphBLAS method, with GxB extensions
GrB_Index Tp_size = (n+1) * sizeof (GrB_Index) ;
GrB_Index Tj_size = nvals * sizeof (GrB_Index) ;
GrB_Index Tx_size = sizeof (bool) ;
Tp = LAGraph_Malloc (n+1, sizeof (GrB_Index)) ;
Tj = LAGraph_Malloc (nvals, sizeof (GrB_Index)) ;
Tx = LAGraph_Calloc (1, sizeof (bool)) ;
range = LAGraph_Malloc (nthreads + 1, sizeof (int64_t)) ;
count = LAGraph_Calloc (nthreads + 1, sizeof (GrB_Index)) ;
Tp = (GrB_Index *) LAGraph_Malloc (n+1, sizeof (GrB_Index)) ;
Tj = (GrB_Index *) LAGraph_Malloc (nvals, sizeof (GrB_Index)) ;
Tx = (bool *) LAGraph_Calloc (1, sizeof (bool)) ;
range = (int64_t *) LAGraph_Malloc (nthreads + 1, sizeof (int64_t)) ;
count = (GrB_Index *) LAGraph_Calloc (nthreads + 1, sizeof (GrB_Index));
LG_CHECK (Tp == NULL || Tj == NULL || Tx == NULL || range == NULL
|| count == NULL, GrB_OUT_OF_MEMORY, "out of memory") ;

Expand Down Expand Up @@ -356,32 +386,31 @@ int LG_CC_FastSV6 // SuiteSparse:GraphBLAS method, with GxB extensions
{
// hooking & shortcutting
// mngp = min (mngp, A*gp) using the MIN_SECOND semiring
GrB_TRY (GrB_mxv (mngp, NULL, GrB_MIN_UINT64,
GrB_MIN_SECOND_SEMIRING_UINT64, T, gp, NULL)) ;
GrB_TRY (GrB_mxv (mngp, NULL, min, min_2nd, T, gp, NULL)) ;

// parent = min (parent, C*mngp) where C is C(i,j) = true if i=Px(j)
GrB_TRY (Reduce_assign (parent, mngp, C, &Cp, &Px, &Cx, msg)) ;
// parent = min (parent, C*mngp), where C(i,j) present if i=Px(j)
GrB_TRY (Reduce_assign (parent, min, min_2nd, mngp, C, &Cp, &Px,
&Cx, msg)) ;

// parent = min (parent, mngp, gp)
GrB_TRY (GrB_eWiseAdd (parent, NULL, GrB_MIN_UINT64, GrB_MIN_UINT64,
mngp, gp, NULL)) ;
GrB_TRY (GrB_eWiseAdd (parent, NULL, min, min, mngp, gp, NULL)) ;

// calculate grandparent: gp_new = parent (parent)
// if parent is uint32, GraphBLAS typecasts to uint64 for Px
GrB_TRY (GrB_Vector_extractTuples (NULL, Px, &n, parent)) ;
GrB_TRY (GrB_extract (gp_new, NULL, NULL, parent, Px, n, NULL)) ;

// terminate if gp and gp_new are the same
GrB_TRY (GrB_eWiseMult (mod, NULL, NULL, GrB_NE_UINT64, gp_new,
gp, NULL)) ;
GrB_TRY (GrB_reduce (&changing, NULL, GrB_LOR_MONOID_BOOL, mod,
GrB_TRY (GrB_eWiseMult (t, NULL, NULL, ne, gp_new, gp, NULL)) ;
GrB_TRY (GrB_reduce (&changing, NULL, GrB_LOR_MONOID_BOOL, t,
NULL)) ;

// swap gp and gp_new
GrB_Vector t = gp ; gp = gp_new ; gp_new = t ;
}

//----------------------------------------------------------------------
// use samping to estimate the largest connected component in T
// use sampling to estimate the largest connected component in T
//----------------------------------------------------------------------

// hash table size must be a power of 2
Expand All @@ -392,26 +421,26 @@ int LG_CC_FastSV6 // SuiteSparse:GraphBLAS method, with GxB extensions
#define NEXT(x) ((x + 23) & (HASH_SIZE-1))

// allocate and initialize the hash table
ht_key = LAGraph_Malloc (HASH_SIZE, sizeof (int64_t)) ;
ht_count = LAGraph_Calloc (HASH_SIZE, sizeof (int64_t)) ;
ht_key = (GrB_Index *) LAGraph_Malloc (HASH_SIZE, sizeof (GrB_Index)) ;
ht_count = (int *) LAGraph_Calloc (HASH_SIZE, sizeof (int)) ;
LG_CHECK (ht_key == NULL || ht_count == NULL, GrB_OUT_OF_MEMORY,
"out of memory") ;
for (int k = 0 ; k < HASH_SIZE ; k++)
{
ht_key [k] = -1 ;
ht_key [k] = UINT64_MAX ;
}

// hash the samples and find the most frequent entry
uint64_t seed = n ; // random number seed
int64_t key = -1 ; // most frequent entry
int64_t max_count = 0 ; // frequency of most frequent entry
int max_count = 0 ; // frequency of most frequent entry
for (int64_t k = 0 ; k < HASH_SAMPLES ; k++)
{
// select an entry from Px at random
int64_t x = Px [LAGraph_Random60 (&seed) % n] ;
GrB_Index x = Px [LAGraph_Random60 (&seed) % n] ;
// find x in the hash table
int64_t h = HASH (x) ;
while (ht_key [h] != -1 && ht_key [h] != x)
GrB_Index h = HASH (x) ;
while (ht_key [h] != UINT64_MAX && ht_key [h] != x)
{
h = NEXT (h) ;
}
Expand All @@ -430,7 +459,7 @@ int LG_CC_FastSV6 // SuiteSparse:GraphBLAS method, with GxB extensions
// compact the largest connected component in T
//----------------------------------------------------------------------

// TODO: replace this with GrB_select and GrB_assign
// TODO: replace this with GxB_extract with GrB_Vector index arrays.

// unpack T to resuse the space (all content is overwritten below)
bool T_jumbled, T_iso ;
Expand Down Expand Up @@ -470,14 +499,14 @@ int LG_CC_FastSV6 // SuiteSparse:GraphBLAS method, with GxB extensions
// key == Px [j]. One of these j's can then be replaced
// with the key. If node i is not adjacent to any node in
// the largest component, then there is no space in T(i,:)
// and no new edge to the largest compenent is added.
// and no new edge to the largest component is added.
if (p - Tp [i] < Sp [i+1] - Sp [i])
{
Tj [p++] = key ;
}
}
}
// count the number of entries inserted into T by this thread?
// count the number of entries inserted into T by this thread
count [tid] = p - Tp [range [tid]] ;
}

Expand Down Expand Up @@ -537,24 +566,23 @@ int LG_CC_FastSV6 // SuiteSparse:GraphBLAS method, with GxB extensions
{
// hooking & shortcutting
// mngp = min (mngp, A*gp) using the MIN_SECOND semiring
GrB_TRY (GrB_mxv (mngp, NULL, GrB_MIN_UINT64,
GrB_MIN_SECOND_SEMIRING_UINT64, A, gp, NULL)) ;
GrB_TRY (GrB_mxv (mngp, NULL, min, min_2nd, A, gp, NULL)) ;

// parent = min (parent, C*mngp) where C is C(i,j) = true if i=Px(j)
GrB_TRY (Reduce_assign (parent, mngp, C, &Cp, &Px, &Cx, msg)) ;
// parent = min (parent, C*mngp) where C(i,j) is present if i=Px(j)
GrB_TRY (Reduce_assign (parent, min, min_2nd, mngp, C, &Cp, &Px, &Cx,
msg)) ;

// parent = min (parent, mngp, gp)
GrB_TRY (GrB_eWiseAdd (parent, NULL, GrB_MIN_UINT64, GrB_MIN_UINT64,
mngp, gp, NULL)) ;
GrB_TRY (GrB_eWiseAdd (parent, NULL, min, min, mngp, gp, NULL)) ;

// calculate grandparent: gp_new = parent (parent)
// if parent is uint32, GraphBLAS typecasts to uint64 for Px.
GrB_TRY (GrB_Vector_extractTuples (NULL, Px, &n, parent)) ;
GrB_TRY (GrB_extract (gp_new, NULL, NULL, parent, Px, n, NULL)) ;

// terminate if gp and gp_new are the same
GrB_TRY (GrB_eWiseMult (mod, NULL, NULL, GrB_NE_UINT64, gp_new, gp,
NULL)) ;
GrB_TRY (GrB_reduce (&changing, NULL, GrB_LOR_MONOID_BOOL, mod, NULL)) ;
GrB_TRY (GrB_eWiseMult (t, NULL, NULL, ne, gp_new, gp, NULL)) ;
GrB_TRY (GrB_reduce (&changing, NULL, GrB_LOR_MONOID_BOOL, t, NULL)) ;

// swap gp and gp_new
GrB_Vector t = gp ; gp = gp_new ; gp_new = t ;
Expand Down
Loading

0 comments on commit 0b0cb0a

Please sign in to comment.