LCOV - code coverage report
Current view: top level - experimental/algorithm - LAGraph_BF_full1.c (source / functions) Hit Total Coverage
Test: LAGraph code coverage report. Commit id: 7b9d2ee. Current time (UTC): 2025-06-03T21:57:17Z Lines: 111 111 100.0 %
Date: 2025-06-03 22:02:40 Functions: 5 5 100.0 %

          Line data    Source code
       1             : //------------------------------------------------------------------------------
       2             : // LAGraph_BF_full1.c: Bellman-Ford single-source shortest paths, returns tree,
       3             : // while diagonal of input matrix A needs not to be explicit 0
       4             : //------------------------------------------------------------------------------
       5             : 
       6             : // LAGraph, (c) 2019-2022 by The LAGraph Contributors, All Rights Reserved.
       7             : // SPDX-License-Identifier: BSD-2-Clause
       8             : //
       9             : // For additional details (including references to third party source code and
      10             : // other files) see the LICENSE file or contact permission@sei.cmu.edu. See
      11             : // Contributors.txt for a full list of contributors. Created, in part, with
      12             : // funding and support from the U.S. Government (see Acknowledgments.txt file).
      13             : // DM22-0790
      14             : 
      15             : // Contributed by Jinhao Chen and Timothy A. Davis, Texas A&M University
      16             : 
      17             : //------------------------------------------------------------------------------
      18             : 
      19             : // LAGraph_BF_full1: Bellman-Ford single source shortest paths, returning both
      20             : // the path lengths and the shortest-path tree.
      21             : 
      22             : // LAGraph_BF_full1 performs a Bellman-Ford to find out shortest path, parent
      23             : // nodes along the path and the hops (number of edges) in the path from given
      24             : // source vertex s in the range of [0, n) on graph given as matrix A with size
      25             : // n*n. The sparse matrix A has entry A(i, j) if there is an edge from vertex i
      26             : // to vertex j with weight w, then A(i, j) = w.
      27             : 
      28             : // LAGraph_BF_full1 returns GrB_SUCCESS if it succeeds.  In this case, there
      29             : // are no negative-weight cycles in the graph, and d, pi, and h are returned.
      30             : // The vector d has d(k) as the shortest distance from s to k. pi(k) = p+1,
      31             : // where p is the parent node of k-th node in the shortest path. In particular,
      32             : // pi(s) = 0. h(k) = hop(s, k), the number of edges from s to k in the shortest
      33             : // path.
      34             : 
      35             : // If the graph has a negative-weight cycle, GrB_NO_VALUE is returned, and the
      36             : // GrB_Vectors d(k), pi(k) and h(k)  (i.e., *pd_output, *ppi_output and
      37             : // *ph_output respectively) will be NULL when negative-weight cycle detected.
      38             : 
      39             : // Otherwise, other errors such as GrB_OUT_OF_MEMORY, GrB_INVALID_OBJECT, and
      40             : // so on, can be returned, if these errors are found by the underlying
      41             : // GrB_* functions.
      42             : 
      43             : //------------------------------------------------------------------------------
      44             : 
      45             : #define LG_FREE_WORK                   \
      46             : {                                      \
      47             :     GrB_free(&d);                      \
      48             :     GrB_free(&dmasked);                \
      49             :     GrB_free(&dless);                  \
      50             :     GrB_free(&Atmp);                   \
      51             :     GrB_free(&BF_Tuple3);              \
      52             :     GrB_free(&BF_lMIN_Tuple3);         \
      53             :     GrB_free(&BF_PLUSrhs_Tuple3);      \
      54             :     GrB_free(&BF_LT_Tuple3);           \
      55             :     GrB_free(&BF_lMIN_Tuple3_Monoid);  \
      56             :     GrB_free(&BF_lMIN_PLUSrhs_Tuple3); \
      57             :     GrB_free(&BF_Identity_Tuple3 );    \
      58             :     LAGraph_Free ((void**)&I, NULL);   \
      59             :     LAGraph_Free ((void**)&J, NULL);   \
      60             :     LAGraph_Free ((void**)&w, NULL);   \
      61             :     LAGraph_Free ((void**)&W, NULL);   \
      62             :     LAGraph_Free ((void**)&h, NULL);   \
      63             :     LAGraph_Free ((void**)&pi, NULL);  \
      64             : }
      65             : 
      66             : #define LG_FREE_ALL                    \
      67             : {                                      \
      68             :     LG_FREE_WORK ;                     \
      69             :     GrB_free (pd_output);              \
      70             :     GrB_free (ppi_output);             \
      71             :     GrB_free (ph_output);              \
      72             : }
      73             : 
      74             : #include <LAGraph.h>
      75             : #include <LAGraphX.h>
      76             : #include <LG_internal.h>  // from src/utility
      77             : 
      78             : 
      79             : typedef void (*LAGraph_binary_function) (void *, const void *, const void *) ;
      80             : 
      81             : //------------------------------------------------------------------------------
      82             : // data type for each entry of the adjacent matrix A and "distance" vector d;
      83             : // <INFINITY,INFINITY,INFINITY> corresponds to nonexistence of a path, and
      84             : // the value  <0, 0, NULL> corresponds to a path from a vertex to itself
      85             : //------------------------------------------------------------------------------
      86             : typedef struct
      87             : {
      88             :     double w;    // w  corresponds to a path weight.
      89             :     GrB_Index h; // h  corresponds to a path size or number of hops.
      90             :     GrB_Index pi;// pi corresponds to the penultimate vertex along a path.
      91             :                  // vertex indexed as 1, 2, 3, ... , V, and pi = 0 (as nil)
      92             :                  // for u=v, and pi = UINT64_MAX (as inf) for (u,v) not in E
      93             : }
      94             : BF1_Tuple3_struct;
      95             : 
      96             : //------------------------------------------------------------------------------
      97             : // 2 binary functions, z=f(x,y), where Tuple3xTuple3 -> Tuple3
      98             : //------------------------------------------------------------------------------
      99       13716 : void BF1_lMIN
     100             : (
     101             :     BF1_Tuple3_struct *z,
     102             :     const BF1_Tuple3_struct *x,
     103             :     const BF1_Tuple3_struct *y
     104             : )
     105             : {
     106       13716 :     if (x->w < y->w
     107        4243 :         || (x->w == y->w && x->h < y->h)
     108        4243 :         || (x->w == y->w && x->h == y->h && x->pi < y->pi))
     109             :     {
     110        9568 :         if (z != x) { *z = *x; }
     111             :     }
     112             :     else
     113             :     {
     114        4148 :         *z = *y;
     115             :     }
     116       13716 : }
     117             : 
     118       18388 : void BF1_PLUSrhs
     119             : (
     120             :     BF1_Tuple3_struct *z,
     121             :     const BF1_Tuple3_struct *x,
     122             :     const BF1_Tuple3_struct *y
     123             : )
     124             : {
     125       18388 :     z->w = x->w + y->w;
     126       18388 :     z->h = x->h + y->h;
     127       18388 :     z->pi = (x->pi != UINT64_MAX && y->pi != 0) ?  y->pi : x->pi ;
     128       18388 : }
     129             : 
     130        9645 : void BF1_Identity
     131             : (
     132             :     BF1_Tuple3_struct *z,
     133             :     const BF1_Tuple3_struct *x
     134             : )
     135             : {
     136        9645 :     *z = *x;
     137        9645 : }
     138             : 
     139        4672 : void BF1_LT
     140             : (
     141             :     bool *z,
     142             :     const BF1_Tuple3_struct *x,
     143             :     const BF1_Tuple3_struct *y
     144             : )
     145             : {
     146        9344 :     (*z) = (x->w < y->w
     147         428 :         || (x->w == y->w && x->h < y->h)
     148        5100 :         || (x->w == y->w && x->h == y->h && x->pi < y->pi)) ;
     149        4672 : }
     150             : 
     151             : // Given a n-by-n adjacency matrix A and a source vertex s.
     152             : // If there is no negative-weight cycle reachable from s, return the distances
     153             : // of shortest paths from s and parents along the paths as vector d. Otherwise,
     154             : // returns d=NULL if there is a negtive-weight cycle.
     155             : // pd_output is pointer to a GrB_Vector, where the i-th entry is d(s,i), the
     156             : //   sum of edges length in the shortest path
     157             : // ppi_output is pointer to a GrB_Vector, where the i-th entry is pi(i), the
     158             : //   parent of i-th vertex in the shortest path
     159             : // ph_output is pointer to a GrB_Vector, where the i-th entry is h(s,i), the
     160             : //   number of edges from s to i in the shortest path
     161             : // A has weights on corresponding entries of edges
     162             : // s is given index for source vertex
     163           5 : GrB_Info LAGraph_BF_full1
     164             : (
     165             :     GrB_Vector *pd_output,      //the pointer to the vector of distance
     166             :     GrB_Vector *ppi_output,     //the pointer to the vector of parent
     167             :     GrB_Vector *ph_output,      //the pointer to the vector of hops
     168             :     const GrB_Matrix A,         //matrix for the graph
     169             :     const GrB_Index s           //given index of the source
     170             : )
     171             : {
     172             :     GrB_Info info;
     173           5 :     char *msg = NULL ;
     174             :     // tmp vector to store distance vector after n (i.e., V) loops
     175           5 :     GrB_Vector d = NULL, dmasked = NULL, dless = NULL;
     176           5 :     GrB_Matrix Atmp = NULL;
     177           5 :     GrB_Type BF_Tuple3 = NULL ;
     178             : 
     179           5 :     GrB_BinaryOp BF_lMIN_Tuple3 = NULL ;
     180           5 :     GrB_BinaryOp BF_PLUSrhs_Tuple3 = NULL ;
     181           5 :     GrB_UnaryOp BF_Identity_Tuple3 = NULL ;
     182           5 :     GrB_BinaryOp BF_LT_Tuple3 = NULL ;
     183             : 
     184           5 :     GrB_Monoid BF_lMIN_Tuple3_Monoid = NULL ;
     185           5 :     GrB_Semiring BF_lMIN_PLUSrhs_Tuple3 = NULL ;
     186             : 
     187             :     GrB_Index nrows, ncols, n, nz;  // n = # of row/col, nz = # of nnz in graph
     188           5 :     GrB_Index *I = NULL, *J = NULL; // for col/row indices of entries from A
     189           5 :     GrB_Index *h = NULL, *pi = NULL;
     190           5 :     double *w = NULL;
     191           5 :     BF1_Tuple3_struct *W = NULL;
     192             : 
     193           5 :     LG_ASSERT (A != NULL && pd_output != NULL &&
     194             :         ppi_output != NULL && ph_output != NULL, GrB_NULL_POINTER) ;
     195             : 
     196           5 :     *pd_output  = NULL;
     197           5 :     *ppi_output = NULL;
     198           5 :     *ph_output  = NULL;
     199           5 :     GRB_TRY (GrB_Matrix_nrows (&nrows, A)) ;
     200           5 :     GRB_TRY (GrB_Matrix_ncols (&ncols, A)) ;
     201           5 :     GRB_TRY (GrB_Matrix_nvals (&nz, A));
     202           5 :     LG_ASSERT_MSG (nrows == ncols, -1002, "A must be square") ;
     203           5 :     n = nrows;
     204           5 :     LG_ASSERT_MSG (s < n, GrB_INVALID_INDEX, "invalid source node") ;
     205             : 
     206             :     //--------------------------------------------------------------------------
     207             :     // create all GrB_Type GrB_BinaryOp GrB_Monoid and GrB_Semiring
     208             :     //--------------------------------------------------------------------------
     209             :     // GrB_Type
     210           5 :     GRB_TRY (GrB_Type_new(&BF_Tuple3, sizeof(BF1_Tuple3_struct)));
     211             : 
     212             :     // GrB_BinaryOp
     213           5 :     GRB_TRY (GrB_UnaryOp_new(&BF_Identity_Tuple3,
     214             :         (void*) (&BF1_Identity), BF_Tuple3, BF_Tuple3));
     215           5 :     GRB_TRY (GrB_BinaryOp_new(&BF_LT_Tuple3,
     216             :         (LAGraph_binary_function) (&BF1_LT), GrB_BOOL, BF_Tuple3, BF_Tuple3));
     217           5 :     GRB_TRY (GrB_BinaryOp_new(&BF_lMIN_Tuple3,
     218             :         (LAGraph_binary_function) (&BF1_lMIN), BF_Tuple3, BF_Tuple3, BF_Tuple3));
     219           5 :     GRB_TRY (GrB_BinaryOp_new(&BF_PLUSrhs_Tuple3,
     220             :         (LAGraph_binary_function)(&BF1_PLUSrhs),
     221             :         BF_Tuple3, BF_Tuple3, BF_Tuple3));
     222             : 
     223             :     // GrB_Monoid
     224           5 :     BF1_Tuple3_struct BF_identity = (BF1_Tuple3_struct) { .w = INFINITY,
     225             :         .h = UINT64_MAX, .pi = UINT64_MAX };
     226           5 :     GRB_TRY (GrB_Monoid_new_UDT(&BF_lMIN_Tuple3_Monoid, BF_lMIN_Tuple3,
     227             :         &BF_identity));
     228             : 
     229             :     //GrB_Semiring
     230           5 :     GRB_TRY (GrB_Semiring_new(&BF_lMIN_PLUSrhs_Tuple3,
     231             :         BF_lMIN_Tuple3_Monoid, BF_PLUSrhs_Tuple3));
     232             : 
     233             :     //--------------------------------------------------------------------------
     234             :     // allocate arrays used for tuplets
     235             :     //--------------------------------------------------------------------------
     236             : 
     237           5 :     LAGRAPH_TRY (LAGraph_Malloc ((void **) &I, nz, sizeof(GrB_Index), msg)) ;
     238           5 :     LAGRAPH_TRY (LAGraph_Malloc ((void **) &J, nz, sizeof(GrB_Index), msg)) ;
     239           5 :     LAGRAPH_TRY (LAGraph_Malloc ((void **) &w, nz, sizeof(double), msg)) ;
     240           5 :     LAGRAPH_TRY (LAGraph_Malloc ((void **) &W, nz, sizeof(BF1_Tuple3_struct),
     241             :         msg)) ;
     242             : 
     243             :     //--------------------------------------------------------------------------
     244             :     // create matrix Atmp based on A, while its entries become BF_Tuple3 type
     245             :     //--------------------------------------------------------------------------
     246             : 
     247           5 :     GRB_TRY (GrB_Matrix_extractTuples_FP64(I, J, w, &nz, A));
     248             :     int nthreads, nthreads_outer, nthreads_inner ;
     249           5 :     LG_TRY (LAGraph_GetNumThreads (&nthreads_outer, &nthreads_inner, msg)) ;
     250           5 :     nthreads = nthreads_outer * nthreads_inner ;
     251             :     int64_t k;
     252             :     #pragma omp parallel for num_threads(nthreads) schedule(static)
     253         773 :     for (k = 0; k < nz; k++)
     254             :     {
     255         768 :         W[k] = (BF1_Tuple3_struct) { .w = w[k], .h = 1, .pi = I[k] + 1 };
     256             :     }
     257           5 :     GRB_TRY (GrB_Matrix_new(&Atmp, BF_Tuple3, n, n));
     258           5 :     GRB_TRY (GrB_Matrix_build_UDT(Atmp, I, J, W, nz, BF_lMIN_Tuple3));
     259           5 :     LAGraph_Free ((void**)&I, NULL);
     260           5 :     LAGraph_Free ((void**)&J, NULL);
     261           5 :     LAGraph_Free ((void**)&W, NULL);
     262           5 :     LAGraph_Free ((void**)&w, NULL);
     263             : 
     264             :     //--------------------------------------------------------------------------
     265             :     // create and initialize "distance" vector d, dmasked and dless
     266             :     //--------------------------------------------------------------------------
     267           5 :     GRB_TRY (GrB_Vector_new(&d, BF_Tuple3, n));
     268             :     // make d dense
     269           5 :     GRB_TRY (GrB_Vector_assign_UDT(d, NULL, NULL, (void*)&BF_identity,
     270             :         GrB_ALL, n, NULL));
     271             :     // initial distance from s to itself
     272           5 :     BF1_Tuple3_struct d0 = (BF1_Tuple3_struct) { .w = 0, .h = 0, .pi = 0 };
     273           5 :     GRB_TRY (GrB_Vector_setElement_UDT(d, &d0, s));
     274             : 
     275             :     // creat dmasked as a sparse vector with only one entry at s
     276           5 :     GRB_TRY (GrB_Vector_new(&dmasked, BF_Tuple3, n));
     277           5 :     GRB_TRY (GrB_Vector_setElement_UDT(dmasked, &d0, s));
     278             : 
     279             :     // create dless
     280           5 :     GRB_TRY (GrB_Vector_new(&dless, GrB_BOOL, n));
     281             : 
     282             :     //--------------------------------------------------------------------------
     283             :     // start the Bellman Ford process
     284             :     //--------------------------------------------------------------------------
     285           5 :     bool any_dless= true;      // if there is any newly found shortest path
     286           5 :     int64_t iter = 0;          // number of iterations
     287             : 
     288             :     // terminate when no new path is found or more than V-1 loops
     289          93 :     while (any_dless && iter < n - 1)
     290             :     {
     291             :         // execute semiring on d and A, and save the result to dtmp
     292          88 :         GRB_TRY (GrB_vxm(dmasked, GrB_NULL, GrB_NULL,
     293             :             BF_lMIN_PLUSrhs_Tuple3, dmasked, Atmp, GrB_NULL));
     294             : 
     295             :         // dless = d .< dtmp
     296             :         //GRB_TRY (GrB_Vector_clear(dless));
     297          88 :         GRB_TRY (GrB_eWiseMult(dless, NULL, NULL, BF_LT_Tuple3, dmasked, d,
     298             :             NULL));
     299             : 
     300             :         // if there is no entry with smaller distance then all shortest paths
     301             :         // are found
     302          88 :         GRB_TRY (GrB_reduce (&any_dless, NULL, GrB_LOR_MONOID_BOOL, dless,
     303             :             NULL)) ;
     304          88 :         if(any_dless)
     305             :         {
     306             :             // update all entries with smaller distances
     307          85 :             GRB_TRY (GrB_apply(d, dless, NULL, BF_Identity_Tuple3, dmasked, NULL));
     308             : 
     309             :             // only use entries that were just updated
     310          85 :             GRB_TRY (GrB_Vector_clear(dmasked));
     311          85 :             GRB_TRY (GrB_apply(dmasked, dless, NULL, BF_Identity_Tuple3, d, NULL));
     312             :             //try:
     313             :             //GRB_TRY (GrB_assign(dmasked, dless, NULL, d, GrB_ALL, n, GrB_DESC_R);
     314             :         }
     315          88 :         iter ++;
     316             :     }
     317             : 
     318             :     // check for negative-weight cycle only when there was a new path in the
     319             :     // last loop, otherwise, there can't be a negative-weight cycle.
     320           5 :     if (any_dless)
     321             :     {
     322             :         // execute semiring again to check for negative-weight cycle
     323           2 :         GRB_TRY (GrB_vxm(dmasked, GrB_NULL, GrB_NULL,
     324             :             BF_lMIN_PLUSrhs_Tuple3, dmasked, Atmp, GrB_NULL));
     325             : 
     326             :         // dless = d .< dtmp
     327             :         //GRB_TRY (GrB_Vector_clear(dless));
     328           2 :         GRB_TRY (GrB_eWiseMult(dless, NULL, NULL, BF_LT_Tuple3, dmasked, d, NULL));
     329             : 
     330             :         // if there is no entry with smaller distance then all shortest paths
     331             :         // are found
     332           2 :         GRB_TRY (GrB_reduce (&any_dless, NULL, GrB_LOR_MONOID_BOOL, dless, NULL)) ;
     333           2 :         if(any_dless)
     334             :         {
     335           2 :             LG_FREE_ALL;
     336           2 :             return (GrB_NO_VALUE) ;
     337             :         }
     338             :     }
     339             : 
     340             :     //--------------------------------------------------------------------------
     341             :     // extract tuple from "distance" vector d and create GrB_Vectors for output
     342             :     //--------------------------------------------------------------------------
     343             : 
     344           3 :     LAGRAPH_TRY (LAGraph_Malloc ((void **) &I, n, sizeof(GrB_Index), msg)) ;
     345           3 :     LAGRAPH_TRY (LAGraph_Malloc ((void **) &W, n, sizeof(BF1_Tuple3_struct),
     346             :         msg)) ;
     347           3 :     LAGRAPH_TRY (LAGraph_Malloc ((void **) &w, n, sizeof(double), msg)) ;
     348           3 :     LAGRAPH_TRY (LAGraph_Malloc ((void **) &h, n, sizeof(GrB_Index), msg)) ;
     349           3 :     LAGRAPH_TRY (LAGraph_Malloc ((void **) &pi, n, sizeof(GrB_Index), msg)) ;
     350             : 
     351           3 :     GRB_TRY (GrB_Vector_extractTuples_UDT (I, (void *) W, &n, d));
     352             : 
     353         111 :     for (k = 0; k < n; k++)
     354             :     {
     355         108 :         w [k] = W[k].w ;
     356         108 :         h [k] = W[k].h ;
     357         108 :         pi[k] = W[k].pi;
     358             :     }
     359           3 :     GRB_TRY (GrB_Vector_new(pd_output,  GrB_FP64,   n));
     360           3 :     GRB_TRY (GrB_Vector_new(ppi_output, GrB_UINT64, n));
     361           3 :     GRB_TRY (GrB_Vector_new(ph_output,  GrB_UINT64, n));
     362           3 :     GRB_TRY (GrB_Vector_build (*pd_output , I, w , n, GrB_MIN_FP64  ));
     363           3 :     GRB_TRY (GrB_Vector_build (*ppi_output, I, pi, n, GrB_MIN_UINT64));
     364           3 :     GRB_TRY (GrB_Vector_build (*ph_output , I, h , n, GrB_MIN_UINT64));
     365           3 :     LG_FREE_WORK;
     366           3 :     return (GrB_SUCCESS) ;
     367             : }

Generated by: LCOV version 1.14