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

          Line data    Source code
       1             : //------------------------------------------------------------------------------
       2             : // LAGr_SingleSourceShortestPath: single-source shortest path
       3             : //------------------------------------------------------------------------------
       4             : 
       5             : // LAGraph, (c) 2019-2022 by The LAGraph Contributors, All Rights Reserved.
       6             : // SPDX-License-Identifier: BSD-2-Clause
       7             : //
       8             : // For additional details (including references to third party source code and
       9             : // other files) see the LICENSE file or contact permission@sei.cmu.edu. See
      10             : // Contributors.txt for a full list of contributors. Created, in part, with
      11             : // funding and support from the U.S. Government (see Acknowledgments.txt file).
      12             : // DM22-0790
      13             : 
      14             : // Contributed by Jinhao Chen, Scott Kolodziej and Tim Davis, Texas A&M
      15             : // University.  Adapted from GraphBLAS Template Library (GBTL) by Scott
      16             : // McMillan and Tze Meng Low.
      17             : 
      18             : //------------------------------------------------------------------------------
      19             : 
      20             : // This is an Advanced algorithm (G->emin is required).
      21             : 
      22             : // Single source shortest path with delta stepping.
      23             : 
      24             : // U. Sridhar, M. Blanco, R. Mayuranath, D. G. Spampinato, T. M. Low, and
      25             : // S. McMillan, "Delta-Stepping SSSP: From Vertices and Edges to GraphBLAS
      26             : // Implementations," in 2019 IEEE International Parallel and Distributed
      27             : // Processing Symposium Workshops (IPDPSW), 2019, pp. 241–250.
      28             : // https://ieeexplore.ieee.org/document/8778222/references
      29             : // https://arxiv.org/abs/1911.06895
      30             : 
      31             : // LAGr_SingleSourceShortestPath computes the shortest path lengths from the
      32             : // specified source vertex to all other vertices in the graph.
      33             : 
      34             : // The parent vector is not computed; see LAGraph_BF_* instead.
      35             : 
      36             : // NOTE: this method gets stuck in an infinite loop when there are negative-
      37             : // weight cycles in the graph.
      38             : 
      39             : // FUTURE: a Basic algorithm that picks Delta automatically
      40             : 
      41             : #define LG_FREE_WORK        \
      42             : {                           \
      43             :     GrB_free (&AL) ;        \
      44             :     GrB_free (&AH) ;        \
      45             :     GrB_free (&lBound) ;    \
      46             :     GrB_free (&uBound) ;    \
      47             :     GrB_free (&tmasked) ;   \
      48             :     GrB_free (&tReq) ;      \
      49             :     GrB_free (&tless) ;     \
      50             :     GrB_free (&s) ;         \
      51             :     GrB_free (&reach) ;     \
      52             :     GrB_free (&Empty) ;     \
      53             : }
      54             : 
      55             : #define LG_FREE_ALL         \
      56             : {                           \
      57             :     LG_FREE_WORK ;          \
      58             :     GrB_free (&t) ;         \
      59             : }
      60             : 
      61             : #include "LG_internal.h"
      62             : 
      63             : #define setelement(s, k)                                                      \
      64             : {                                                                             \
      65             :     switch (tcode)                                                            \
      66             :     {                                                                         \
      67             :         default:                                                              \
      68             :         case 0 : GrB_Scalar_setElement_INT32  (s, k * delta_int32 ) ; break ; \
      69             :         case 1 : GrB_Scalar_setElement_INT64  (s, k * delta_int64 ) ; break ; \
      70             :         case 2 : GrB_Scalar_setElement_UINT32 (s, k * delta_uint32) ; break ; \
      71             :         case 3 : GrB_Scalar_setElement_UINT64 (s, k * delta_uint64) ; break ; \
      72             :         case 4 : GrB_Scalar_setElement_FP32   (s, k * delta_fp32  ) ; break ; \
      73             :         case 5 : GrB_Scalar_setElement_FP64   (s, k * delta_fp64  ) ; break ; \
      74             :     }                                                                         \
      75             : }
      76             : 
      77        4437 : int LAGr_SingleSourceShortestPath
      78             : (
      79             :     // output:
      80             :     GrB_Vector *path_length,    // path_length (i) is the length of the shortest
      81             :                                 // path from the source vertex to vertex i
      82             :     // input:
      83             :     const LAGraph_Graph G,      // input graph, not modified
      84             :     GrB_Index source,           // source vertex
      85             :     GrB_Scalar Delta,           // delta value for delta stepping
      86             :     char *msg
      87             : )
      88             : {
      89             : 
      90             :     //--------------------------------------------------------------------------
      91             :     // check inputs
      92             :     //--------------------------------------------------------------------------
      93             : 
      94        4437 :     LG_CLEAR_MSG ;
      95        4437 :     GrB_Scalar lBound = NULL ;  // the threshold for GrB_select
      96        4437 :     GrB_Scalar uBound = NULL ;  // the threshold for GrB_select
      97        4437 :     GrB_Matrix AL = NULL ;      // graph containing the light weight edges
      98        4437 :     GrB_Matrix AH = NULL ;      // graph containing the heavy weight edges
      99        4437 :     GrB_Vector t = NULL ;       // tentative shortest path length
     100        4437 :     GrB_Vector tmasked = NULL ;
     101        4437 :     GrB_Vector tReq = NULL ;
     102        4437 :     GrB_Vector tless = NULL ;
     103        4437 :     GrB_Vector s = NULL ;
     104        4437 :     GrB_Vector reach = NULL ;
     105        4437 :     GrB_Vector Empty = NULL ;
     106             : 
     107        4437 :     LG_TRY (LAGraph_CheckGraph (G, msg)) ;
     108        4437 :     LG_ASSERT (path_length != NULL && Delta != NULL, GrB_NULL_POINTER) ;
     109        4437 :     (*path_length) = NULL ;
     110             : 
     111             :     GrB_Index nvals ;
     112        4437 :     LG_TRY (GrB_Scalar_nvals (&nvals, Delta)) ;
     113        4437 :     LG_ASSERT_MSG (nvals == 1, GrB_EMPTY_OBJECT, "Delta is missing") ;
     114             : 
     115        4436 :     GrB_Matrix A = G->A ;
     116             :     GrB_Index n ;
     117        4436 :     GRB_TRY (GrB_Matrix_nrows (&n, A)) ;
     118        4436 :     LG_ASSERT_MSG (source < n, GrB_INVALID_INDEX, "invalid source node") ;
     119             : 
     120             :     //--------------------------------------------------------------------------
     121             :     // initializations
     122             :     //--------------------------------------------------------------------------
     123             : 
     124             :     // get the type of the A matrix
     125             :     GrB_Type etype ;
     126             :     char typename [LAGRAPH_MAX_NAME_LEN] ;
     127        4436 :     LG_TRY (LAGraph_Matrix_TypeName (typename, A, msg)) ;
     128        4436 :     LG_TRY (LAGraph_TypeFromName (&etype, typename, msg)) ;
     129             : 
     130        4436 :     GRB_TRY (GrB_Scalar_new (&lBound, etype)) ;
     131        4408 :     GRB_TRY (GrB_Scalar_new (&uBound, etype)) ;
     132        4380 :     GRB_TRY (GrB_Vector_new (&t, etype, n)) ;
     133        4352 :     GRB_TRY (GrB_Vector_new (&tmasked, etype, n)) ;
     134        4324 :     GRB_TRY (GrB_Vector_new (&tReq, etype, n)) ;
     135        4296 :     GRB_TRY (GrB_Vector_new (&Empty, GrB_BOOL, n)) ;
     136        4268 :     GRB_TRY (GrB_Vector_new (&tless, GrB_BOOL, n)) ;
     137        4240 :     GRB_TRY (GrB_Vector_new (&s, GrB_BOOL, n)) ;
     138        4212 :     GRB_TRY (GrB_Vector_new (&reach, GrB_BOOL, n)) ;
     139             : 
     140             :     // optional hints for SuiteSparse:GraphBLAS
     141        4184 :     GRB_TRY (LG_SET_FORMAT_HINT (t, LG_BITMAP)) ;
     142        4156 :     GRB_TRY (LG_SET_FORMAT_HINT (tmasked, LG_SPARSE)) ;
     143        4156 :     GRB_TRY (LG_SET_FORMAT_HINT (tReq, LG_SPARSE)) ;
     144        4156 :     GRB_TRY (LG_SET_FORMAT_HINT (tless, LG_SPARSE)) ;
     145        4156 :     GRB_TRY (LG_SET_FORMAT_HINT (s, LG_SPARSE)) ;
     146        4156 :     GRB_TRY (LG_SET_FORMAT_HINT (reach, LG_BITMAP)) ;
     147             : 
     148             :     // select the operators, and set t (:) = infinity
     149             :     GrB_IndexUnaryOp ne, le, ge, lt, gt ;
     150             :     GrB_BinaryOp less_than ;
     151             :     GrB_Semiring min_plus ;
     152             :     int tcode ;
     153             :     int32_t  delta_int32  ;
     154             :     int64_t  delta_int64  ;
     155             :     uint32_t delta_uint32 ;
     156             :     uint64_t delta_uint64 ;
     157             :     float    delta_fp32   ;
     158             :     double   delta_fp64   ;
     159             : 
     160        4128 :     bool negative_edge_weights = true ;
     161             : 
     162        4128 :     if (etype == GrB_INT32)
     163             :     {
     164        3901 :         GRB_TRY (GrB_Scalar_extractElement (&delta_int32, Delta)) ;
     165        3901 :         GRB_TRY (GrB_assign (t, NULL, NULL, (int32_t) INT32_MAX,
     166             :             GrB_ALL, n, NULL)) ;
     167        3873 :         le = GrB_VALUELE_INT32 ;
     168        3873 :         ge = GrB_VALUEGE_INT32 ;
     169        3873 :         lt = GrB_VALUELT_INT32 ;
     170        3873 :         gt = GrB_VALUEGT_INT32 ;
     171        3873 :         less_than = GrB_LT_INT32 ;
     172        3873 :         min_plus = GrB_MIN_PLUS_SEMIRING_INT32 ;
     173        3873 :         tcode = 0 ;
     174             :     }
     175         227 :     else if (etype == GrB_INT64)
     176             :     {
     177          48 :         GRB_TRY (GrB_Scalar_extractElement (&delta_int64, Delta)) ;
     178          48 :         GRB_TRY (GrB_assign (t, NULL, NULL, (int64_t) INT64_MAX,
     179             :             GrB_ALL, n, NULL)) ;
     180          48 :         le = GrB_VALUELE_INT64 ;
     181          48 :         ge = GrB_VALUEGE_INT64 ;
     182          48 :         lt = GrB_VALUELT_INT64 ;
     183          48 :         gt = GrB_VALUEGT_INT64 ;
     184          48 :         less_than = GrB_LT_INT64 ;
     185          48 :         min_plus = GrB_MIN_PLUS_SEMIRING_INT64 ;
     186          48 :         tcode = 1 ;
     187             :     }
     188         179 :     else if (etype == GrB_UINT32)
     189             :     {
     190          12 :         GRB_TRY (GrB_Scalar_extractElement (&delta_uint32, Delta)) ;
     191          12 :         GRB_TRY (GrB_assign (t, NULL, NULL, (uint32_t) UINT32_MAX,
     192             :             GrB_ALL, n, NULL)) ;
     193          12 :         le = GrB_VALUELE_UINT32 ;
     194          12 :         ge = GrB_VALUEGE_UINT32 ;
     195          12 :         lt = GrB_VALUELT_UINT32 ;
     196          12 :         gt = GrB_VALUEGT_UINT32 ;
     197          12 :         less_than = GrB_LT_UINT32 ;
     198          12 :         min_plus = GrB_MIN_PLUS_SEMIRING_UINT32 ;
     199          12 :         tcode = 2 ;
     200          12 :         negative_edge_weights = false ;
     201             :     }
     202         167 :     else if (etype == GrB_UINT64)
     203             :     {
     204          12 :         GRB_TRY (GrB_Scalar_extractElement (&delta_uint64, Delta)) ;
     205          12 :         GRB_TRY (GrB_assign (t, NULL, NULL, (uint64_t) UINT64_MAX,
     206             :             GrB_ALL, n, NULL)) ;
     207          12 :         le = GrB_VALUELE_UINT64 ;
     208          12 :         ge = GrB_VALUEGE_UINT64 ;
     209          12 :         lt = GrB_VALUELT_UINT64 ;
     210          12 :         gt = GrB_VALUEGT_UINT64 ;
     211          12 :         less_than = GrB_LT_UINT64 ;
     212          12 :         min_plus = GrB_MIN_PLUS_SEMIRING_UINT64 ;
     213          12 :         tcode = 3 ;
     214          12 :         negative_edge_weights = false ;
     215             :     }
     216         155 :     else if (etype == GrB_FP32)
     217             :     {
     218           9 :         GRB_TRY (GrB_Scalar_extractElement (&delta_fp32, Delta)) ;
     219           9 :         GRB_TRY (GrB_assign (t, NULL, NULL, (float) INFINITY,
     220             :             GrB_ALL, n, NULL)) ;
     221           9 :         le = GrB_VALUELE_FP32 ;
     222           9 :         ge = GrB_VALUEGE_FP32 ;
     223           9 :         lt = GrB_VALUELT_FP32 ;
     224           9 :         gt = GrB_VALUEGT_FP32 ;
     225           9 :         less_than = GrB_LT_FP32 ;
     226           9 :         min_plus = GrB_MIN_PLUS_SEMIRING_FP32 ;
     227           9 :         tcode = 4 ;
     228             :     }
     229         146 :     else if (etype == GrB_FP64)
     230             :     {
     231         145 :         GRB_TRY (GrB_Scalar_extractElement (&delta_fp64, Delta)) ;
     232         145 :         GRB_TRY (GrB_assign (t, NULL, NULL, (double) INFINITY,
     233             :             GrB_ALL, n, NULL)) ;
     234         145 :         le = GrB_VALUELE_FP64 ;
     235         145 :         ge = GrB_VALUEGE_FP64 ;
     236         145 :         lt = GrB_VALUELT_FP64 ;
     237         145 :         gt = GrB_VALUEGT_FP64 ;
     238         145 :         less_than = GrB_LT_FP64 ;
     239         145 :         min_plus = GrB_MIN_PLUS_SEMIRING_FP64 ;
     240         145 :         tcode = 5 ;
     241             :     }
     242             :     else
     243             :     {
     244           1 :         LG_ASSERT_MSG (false, GrB_NOT_IMPLEMENTED, "type not supported") ;
     245             :     }
     246             : 
     247             :     // check if the graph might have negative edge weights
     248        4099 :     if (negative_edge_weights)
     249             :     {
     250        4075 :         double emin = -1 ;
     251        4075 :         if (G->emin != NULL &&
     252         452 :             (G->emin_state == LAGraph_VALUE ||
     253         238 :              G->emin_state == LAGraph_BOUND))
     254             :         {
     255         214 :             GRB_TRY (GrB_Scalar_extractElement_FP64 (&emin, G->emin)) ;
     256             :         }
     257             : //      else
     258             : //      {
     259             : //          // a future Basic algorithm should compute G->emin and perhaps
     260             : //          // G->emax, then compute Delta automatically.
     261             : //      }
     262        4075 :         negative_edge_weights = (emin < 0) ;
     263             :     }
     264             : 
     265             :     // t (src) = 0
     266        4099 :     GRB_TRY (GrB_Vector_setElement (t, 0, source)) ;
     267             : 
     268             :     // reach (src) = true
     269        4085 :     GRB_TRY (GrB_Vector_setElement (reach, true, source)) ;
     270             : 
     271             :     // s (src) = true
     272        4071 :     GRB_TRY (GrB_Vector_setElement (s, true, source)) ;
     273             : 
     274             :     // AL = A .* (A <= Delta)
     275        4029 :     GRB_TRY (GrB_Matrix_new (&AL, etype, n, n)) ;
     276        3987 :     GRB_TRY (GrB_select (AL, NULL, NULL, le, A, Delta, NULL)) ;
     277        3908 :     GRB_TRY (GrB_wait (AL, GrB_MATERIALIZE)) ;
     278             : 
     279             :     // FUTURE: costly for some problems, taking up to 50% of the total time:
     280             :     // AH = A .* (A > Delta)
     281        3908 :     GRB_TRY (GrB_Matrix_new (&AH, etype, n, n)) ;
     282        3866 :     GRB_TRY (GrB_select (AH, NULL, NULL, gt, A, Delta, NULL)) ;
     283        3762 :     GRB_TRY (GrB_wait (AH, GrB_MATERIALIZE)) ;
     284             : 
     285             :     //--------------------------------------------------------------------------
     286             :     // while (t >= step*Delta) not empty
     287             :     //--------------------------------------------------------------------------
     288             : 
     289        3762 :     for (int64_t step = 0 ; ; step++)
     290        9807 :     {
     291             : 
     292             :         //----------------------------------------------------------------------
     293             :         // tmasked = all entries in t<reach> that are less than (step+1)*Delta
     294             :         //----------------------------------------------------------------------
     295             : 
     296       13569 :         setelement (uBound, (step+1)) ;        // uBound = (step+1) * Delta
     297       16682 :         GRB_TRY (GrB_Vector_clear (tmasked)) ;
     298             : 
     299             :         // tmasked<reach> = t
     300             :         // FUTURE: this is costly, typically using Method 06s in SuiteSparse,
     301             :         // which is a very general-purpose one.  Write a specialized kernel to
     302             :         // exploit the fact that reach and t are bitmap and tmasked starts
     303             :         // empty, or fuse this assignment with the GrB_select below.
     304       13439 :         GRB_TRY (GrB_assign (tmasked, reach, NULL, t, GrB_ALL, n, NULL)) ;
     305             :         // tmasked = select (tmasked < (step+1)*Delta)
     306       13171 :         GRB_TRY (GrB_select (tmasked, NULL, NULL, lt, tmasked, uBound, NULL)) ;
     307             :         // --- alternative:
     308             :         // FUTURE this is slower than the above but should be much faster.
     309             :         // GrB_select is computing a bitmap result then converting it to
     310             :         // sparse.  t and reach are both bitmap and tmasked finally sparse.
     311             :         // tmasked<reach> = select (t < (step+1)*Delta)
     312             :         // GRB_TRY (GrB_select (tmasked, reach, NULL, lt, t, uBound, NULL)) ;
     313             : 
     314             :         GrB_Index tmasked_nvals ;
     315       12939 :         GRB_TRY (GrB_Vector_nvals (&tmasked_nvals, tmasked)) ;
     316             : 
     317             :         //----------------------------------------------------------------------
     318             :         // continue while the current bucket (tmasked) is not empty
     319             :         //----------------------------------------------------------------------
     320             : 
     321       26983 :         while (tmasked_nvals > 0)
     322             :         {
     323             :             // tReq = AL'*tmasked using the min_plus semiring
     324       21575 :             GRB_TRY (GrB_vxm (tReq, NULL, NULL, min_plus, tmasked, AL, NULL)) ;
     325             : 
     326             :             // s<struct(tmasked)> = true
     327       19572 :             GRB_TRY (GrB_assign (s, tmasked, NULL, (bool) true, GrB_ALL, n,
     328             :                 GrB_DESC_S)) ;
     329             : 
     330             :             // if nvals (tReq) is 0, no need to continue the rest of this loop
     331             :             GrB_Index tReq_nvals ;
     332       19365 :             GRB_TRY (GrB_Vector_nvals (&tReq_nvals, tReq)) ;
     333       21590 :             if (tReq_nvals == 0) break ;
     334             : 
     335             :             // tless = (tReq .< t) using set intersection
     336       17735 :             GRB_TRY (GrB_eWiseMult (tless, NULL, NULL, less_than, tReq, t,
     337             :                 NULL)) ;
     338             : 
     339             :             // remove explicit zeros from tless so it can be used as a
     340             :             // structural mask
     341             :             GrB_Index tless_nvals ;
     342       17522 :             GRB_TRY (GrB_select (tless, NULL, NULL, GrB_VALUENE_BOOL,
     343             :                 tless, false, NULL)) ;
     344       17177 :             GRB_TRY (GrB_Vector_nvals (&tless_nvals, tless)) ;
     345       17177 :             if (tless_nvals == 0) break ;
     346             : 
     347             :             // update reachable node list/mask
     348             :             // reach<struct(tless)> = true
     349       14952 :             GRB_TRY (GrB_assign (reach, tless, NULL, (bool) true, GrB_ALL, n,
     350             :                 GrB_DESC_S)) ;
     351             : 
     352             :             // tmasked<struct(tless)> = select (tReq < (step+1)*Delta)
     353       14952 :             GRB_TRY (GrB_Vector_clear (tmasked)) ;
     354       14904 :             GRB_TRY (GrB_select (tmasked, tless, NULL, lt, tReq, uBound,
     355             :                 GrB_DESC_S)) ;
     356             : 
     357             :             // For general graph with some negative weights:
     358       14404 :             if (negative_edge_weights)
     359             :             {
     360             :                 // If all entries of the graph are known to be positive, and
     361             :                 // the entries of tmasked are at least step*Delta, tReq =
     362             :                 // tmasked min.+ AL must be >= step*Delta.  Therefore, there is
     363             :                 // no need to perform this GrB_select with ge to find tmasked
     364             :                 // >= step*Delta from tReq.
     365       11533 :                 setelement (lBound, (step)) ;  // lBound = step*Delta
     366             :                 // tmasked = select entries in tmasked that are >= step*Delta
     367       11533 :                 GRB_TRY (GrB_select (tmasked, NULL, NULL, ge, tmasked, lBound,
     368             :                     NULL)) ;
     369             :             }
     370             : 
     371             :             // t<struct(tless)> = tReq
     372       14044 :             GRB_TRY (GrB_assign (t, tless, NULL, tReq, GrB_ALL, n, GrB_DESC_S));
     373       14044 :             GRB_TRY (GrB_Vector_nvals (&tmasked_nvals, tmasked)) ;
     374             :         }
     375             : 
     376             :         // tmasked<s> = t
     377       10936 :         GRB_TRY (GrB_Vector_clear (tmasked)) ;
     378       10904 :         GRB_TRY (GrB_assign (tmasked, s, NULL, t, GrB_ALL, n, GrB_DESC_S)) ;
     379             : 
     380             :         // tReq = AH'*tmasked using the min_plus semiring
     381       10670 :         GRB_TRY (GrB_vxm (tReq, NULL, NULL, min_plus, tmasked, AH, NULL)) ;
     382             : 
     383             :         // tless = (tReq .< t) using set intersection
     384       10446 :         GRB_TRY (GrB_eWiseMult (tless, NULL, NULL, less_than, tReq, t, NULL)) ;
     385             : 
     386             :         // t<tless> = tReq, which computes t = min (t, tReq)
     387       10350 :         GRB_TRY (GrB_assign (t, tless, NULL, tReq, GrB_ALL, n, NULL)) ;
     388             : 
     389             :         //----------------------------------------------------------------------
     390             :         // find out how many left to be computed
     391             :         //----------------------------------------------------------------------
     392             : 
     393             :         // update reachable node list
     394             :         // reach<tless> = true
     395       10350 :         GRB_TRY (GrB_assign (reach, tless, NULL, (bool) true, GrB_ALL, n,
     396             :             NULL)) ;
     397             : 
     398             :         // remove previous buckets
     399             :         // reach<struct(s)> = Empty
     400       10350 :         GRB_TRY (GrB_assign (reach, s, NULL, Empty, GrB_ALL, n, GrB_DESC_S)) ;
     401             :         GrB_Index nreach ;
     402       10344 :         GRB_TRY (GrB_Vector_nvals (&nreach, reach)) ;
     403       10344 :         if (nreach == 0) break ;
     404             : 
     405        9825 :         GRB_TRY (GrB_Vector_clear (s)) ; // clear s for the next iteration
     406             :     }
     407             : 
     408             :     //--------------------------------------------------------------------------
     409             :     // free workspace and return result
     410             :     //--------------------------------------------------------------------------
     411             : 
     412         519 :     (*path_length) = t ;
     413         519 :     LG_FREE_WORK ;
     414         519 :     return (GrB_SUCCESS) ;
     415             : }

Generated by: LCOV version 1.14