pathrex-sys 0.1.0

//------------------------------------------------------------------------------
// LAGraph_scc.c
//------------------------------------------------------------------------------

// LAGraph, (c) 2019-2022 by The LAGraph Contributors, All Rights Reserved.
// SPDX-License-Identifier: BSD-2-Clause
//
// For additional details (including references to third party source code and
// other files) see the LICENSE file or contact permission@sei.cmu.edu. See
// Contributors.txt for a full list of contributors. Created, in part, with
// funding and support from the U.S. Government (see Acknowledgments.txt file).
// DM22-0790

// Contributed by Yongzhe Zhang (zyz915@gmail.com)

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

// TODO: not ready for src; uses global variables

/**
 * Code is based on the Min-Label algorithm described in the following paper:
 * D. Yan, J. Cheng, K. Xin, Y. Lu, W. Ng, Y. Bu, "Pregel Algorithms for Graph
 * Connectivity Problems with Performance Guarantees"
 * Proc. VLDB Endow. 7, 14 (October 2014), 1821–1832.
 * DOI: https://doi.org/10.14778/2733085.2733089
 **/

#define LG_FREE_ALL ;

#include "LG_internal.h"
#include <LAGraph.h>
#include <LAGraphX.h>

#if LAGRAPH_SUITESPARSE

//****************************************************************************
// global C arrays used in SelectOp
GrB_Index *I = NULL, *V = NULL, *F = NULL, *B = NULL, *M = NULL;

// edge_removal:
//  - remove the edges connected to newly identified SCCs (vertices u with M[u]==1)
//  - remove the edges (u, v) where u and v can never be in the same SCC.
//
// Here's a brief explanation of the second case. After the forward and backward
// propagation, each vertex u has two labels
//  - F[u]: the smallest vertex that can reach u
//  - B[u]: the smallest vertex that is reachable from u
// If two vertices u and v are in the same SCC, then F[u]==F[v] and B[u]==B[v] must
// hold. The converse is not true unless F[u]==B[u]. However, we can safely remove
// an edge (u, v) if either F[u]!=F[v] or B[u]!=B[v] holds, which can accelerate
// the SCC computation in the future rounds.

void edge_removal (bool *z, const void *x, GrB_Index i, GrB_Index j, const void *thunk) ;
void edge_removal (bool *z, const void *x, GrB_Index i, GrB_Index j, const void *thunk)
{
    (*z) = (!M[i] && !M[j] && F[i] == F[j] && B[i] == B[j]) ;
}

//****************************************************************************
// trim_one: remove the edges connected to trivial SCCs
//  - A vertex is a trivial SCC if it has no incoming or outgoing edges.
//  - M[i] = i   | if vertex i is a trivial SCC
//    M[i] = n   | otherwise

void trim_one (bool *z, const void *x, GrB_Index i, GrB_Index j, const void *thunk) ;
void trim_one (bool *z, const void *x, GrB_Index i, GrB_Index j, const void *thunk)
{
    (*z) = (M[i] == M[j]) ;
}

//****************************************************************************
// label propagation
//  - label  : (input/output) labels
//  - mask   : (input) mask
//  - A      : (input) original matrix
//  - AT     : (input) transposed matrix
//  - n      : (input) number of vertices

#undef  LG_FREE_ALL
#define LG_FREE_ALL    \
    GrB_free (&s) ;    \
    GrB_free (&t) ;

static GrB_Info propagate (GrB_Vector label, GrB_Vector mask,
        GrB_Matrix A, GrB_Matrix AT, GrB_Index n, char *msg)
{
    GrB_Info info;
    // semirings

    GrB_Vector s = NULL, t = NULL;
    GRB_TRY (GrB_Vector_new (&s, GrB_UINT64, n));
    GRB_TRY (GrB_Vector_new (&t, GrB_UINT64, n));
    GRB_TRY (GrB_assign (s, mask, 0, label, GrB_ALL, 0, 0));
    GRB_TRY (GrB_assign (t, 0, 0, label, GrB_ALL, 0, 0));

    GrB_Index active;
    while (true)
    {
        GRB_TRY (GrB_vxm (t, 0, GrB_MIN_UINT64,
                                 GrB_MIN_FIRST_SEMIRING_UINT64, s, A, 0));
        GRB_TRY (GrB_eWiseMult (mask, 0, 0, GxB_ISNE_UINT64, t, label, 0));
        GRB_TRY (GrB_assign (label, mask, 0, t, GrB_ALL, 0, 0));
        GRB_TRY (GrB_reduce (&active, 0, GrB_PLUS_MONOID_UINT64, mask, 0));
        if (active == 0) break;
        GRB_TRY (GrB_Vector_clear (s));
        GRB_TRY (GrB_assign (s, mask, 0, label, GrB_ALL, 0, 0));
    }

    LG_FREE_ALL ;
    return GrB_SUCCESS;
}

//****************************************************************************

#undef  LG_FREE_ALL
#define LG_FREE_ALL                         \
    LAGraph_Free ((void **) &I, msg);       \
    LAGraph_Free ((void **) &V, msg);       \
    LAGraph_Free ((void **) &F, msg);       \
    LAGraph_Free ((void **) &B, msg);       \
    LAGraph_Free ((void **) &M, msg);       \
    GrB_free (&ind);                        \
    GrB_free (&inf);                        \
    GrB_free (&f);                          \
    GrB_free (&b);                          \
    GrB_free (&mask);                       \
    GrB_free (&FW);                         \
    GrB_free (&BW);                         \
    GrB_free (&sel1);                       \
    GrB_free (&sel2);                       \
    GrB_free (&scc);

#endif

//****************************************************************************
int LAGraph_scc
(
    GrB_Vector *result,     // output: array of component identifiers
    GrB_Matrix A,           // input matrix
    char *msg
)
{
#if LAGRAPH_SUITESPARSE

    LG_CLEAR_MSG ;

    GrB_Info info;
    GrB_Vector scc = NULL ;
    GrB_Vector ind = NULL ;
    GrB_Vector inf = NULL ;
    GrB_Vector f = NULL, b = NULL, mask = NULL ;
    GrB_IndexUnaryOp sel1 = NULL, sel2 = NULL ;
    GrB_Monoid Add = NULL ;
    GrB_Matrix FW = NULL, BW = NULL;

    if (result == NULL || A == NULL) return (GrB_NULL_POINTER) ;

    GrB_Index n, ncols, nvals;
    GRB_TRY (GrB_Matrix_nrows (&n, A));
    GRB_TRY (GrB_Matrix_ncols (&ncols, A));
    if (n != ncols) return (GrB_DIMENSION_MISMATCH) ;

    // store the graph in both directions (forward / backward)
    GRB_TRY (GrB_Matrix_new (&FW, GrB_BOOL, n, n));
    GRB_TRY (GrB_Matrix_new (&BW, GrB_BOOL, n, n));
    GRB_TRY (GrB_transpose (FW, 0, 0, A, GrB_DESC_T0)); // FW = A
    GRB_TRY (GrB_transpose (BW, 0, 0, A, 0));     // BW = A'

    // check format
    int32_t A_format, AT_format;
    GRB_TRY (GrB_get (FW, &A_format , GrB_STORAGE_ORIENTATION_HINT));
    GRB_TRY (GrB_get (BW, &AT_format, GrB_STORAGE_ORIENTATION_HINT));

    bool is_csr = (A_format == GrB_ROWMAJOR && AT_format == GrB_ROWMAJOR);
    if (!is_csr) return (GrB_INVALID_VALUE) ;

    LG_TRY (LAGraph_Malloc ((void **) &I, n, sizeof (GrB_Index), msg)) ;
    LG_TRY (LAGraph_Malloc ((void **) &V, n, sizeof (GrB_Index), msg)) ;
    LG_TRY (LAGraph_Malloc ((void **) &F, n, sizeof (GrB_Index), msg)) ;
    LG_TRY (LAGraph_Malloc ((void **) &B, n, sizeof (GrB_Index), msg)) ;
    LG_TRY (LAGraph_Malloc ((void **) &M, n, sizeof (GrB_Index), msg)) ;

    for (GrB_Index i = 0; i < n; i++)
        I[i] = V[i] = i;

    // scc: the SCC identifier for each vertex
    // scc[u] == n: not assigned yet
    GRB_TRY (GrB_Vector_new (&scc, GrB_UINT64, n));
    // vector of indices: ind[i] == i
    GRB_TRY (GrB_Vector_new (&ind, GrB_UINT64, n));
    GRB_TRY (GrB_Vector_build (ind, I, V, n, GrB_PLUS_UINT64));
    // vector of infinite value: inf[i] == n
    GRB_TRY (GrB_Vector_new (&inf, GrB_UINT64, n));
    GRB_TRY (GrB_assign (inf, 0, 0, n, GrB_ALL, 0, 0));
    // other vectors
    GRB_TRY (GrB_Vector_new (&f, GrB_UINT64, n));
    GRB_TRY (GrB_Vector_new (&b, GrB_UINT64, n));
    GRB_TRY (GrB_Vector_new (&mask, GrB_UINT64, n));
    GRB_TRY (GrB_IndexUnaryOp_new (&sel1, (void *) trim_one, GrB_BOOL, GrB_UINT64, GrB_UINT64));
    GRB_TRY (GrB_IndexUnaryOp_new (&sel2, (void *) edge_removal, GrB_BOOL, GrB_UINT64, GrB_UINT64));

    // remove trivial SCCs
    GRB_TRY (GrB_reduce (f, 0, GrB_PLUS_UINT64, GrB_PLUS_UINT64, FW, 0));
    GRB_TRY (GrB_reduce (b, 0, GrB_PLUS_UINT64, GrB_PLUS_UINT64, BW, 0));
    GRB_TRY (GrB_eWiseMult (mask, 0, GxB_LAND_UINT64, GxB_LAND_UINT64, f, b, 0));
    GRB_TRY (GrB_Vector_nvals (&nvals, mask));

    GRB_TRY (GrB_assign (scc, 0, 0, ind, GrB_ALL, 0, 0));
    GRB_TRY (GrB_assign (scc, mask, 0, n, GrB_ALL, 0, 0));
    GRB_TRY (GrB_Vector_clear (mask));

    if (nvals < n)
    {
        GRB_TRY (GrB_Vector_extractTuples (I, M, &n, scc));
        GRB_TRY (GrB_select (FW, 0, 0, sel1, FW, 0, 0));
        GRB_TRY (GrB_select (BW, 0, 0, sel1, BW, 0, 0));
    }

    GRB_TRY (GrB_Matrix_nvals (&nvals, FW));
    while (nvals > 0)
    {
        GRB_TRY (GrB_eWiseMult (mask, 0, 0, GxB_ISEQ_UINT64, scc, inf, 0));
        GRB_TRY (GrB_assign (f, 0, 0, ind, GrB_ALL, 0, 0));
        LG_TRY (propagate (f, mask, FW, BW, n, msg));

        GRB_TRY (GrB_eWiseMult (mask, 0, 0, GxB_ISEQ_UINT64, f, ind, 0));
        GRB_TRY (GrB_assign (b, 0, 0, inf, GrB_ALL, 0, 0));
        GRB_TRY (GrB_assign (b, mask, 0, ind, GrB_ALL, 0, 0));
        LG_TRY (propagate (b, mask, BW, FW, n, msg));

        GRB_TRY (GrB_eWiseMult (mask, 0, 0, GxB_ISEQ_UINT64, f, b, 0));
        GRB_TRY (GrB_assign (scc, mask, GrB_MIN_UINT64, f, GrB_ALL, 0, 0));

        GRB_TRY (GrB_Vector_extractTuples (I, F, &n, f));
        GRB_TRY (GrB_Vector_extractTuples (I, B, &n, b));
        GRB_TRY (GrB_Vector_extractTuples (I, M, &n, mask));

        GRB_TRY (GrB_select (FW, 0, 0, sel2, FW, 0, 0));
        GRB_TRY (GrB_select (BW, 0, 0, sel2, BW, 0, 0));

        GRB_TRY (GrB_Matrix_nvals (&nvals, FW));
    }
    GRB_TRY (GrB_eWiseMult (mask, 0, 0, GxB_ISEQ_UINT64, scc, inf, 0));
    GRB_TRY (GrB_assign (scc, mask, 0, ind, GrB_ALL, 0, 0));

    GRB_TRY (GrB_eWiseMult (mask, 0, 0, GxB_ISEQ_UINT64, scc, ind, 0));
    GRB_TRY (GrB_reduce (&nvals, 0, GrB_PLUS_MONOID_UINT64, mask, 0));

    *result = scc;
    scc = NULL;

    LG_FREE_ALL ;
    return (GrB_SUCCESS) ;
#else
    return (GrB_NOT_IMPLEMENTED) ;
#endif
}