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 3955 : 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 3955 : LG_CLEAR_MSG ;
95 3955 : GrB_Scalar lBound = NULL ; // the threshold for GrB_select
96 3955 : GrB_Scalar uBound = NULL ; // the threshold for GrB_select
97 3955 : GrB_Matrix AL = NULL ; // graph containing the light weight edges
98 3955 : GrB_Matrix AH = NULL ; // graph containing the heavy weight edges
99 3955 : GrB_Vector t = NULL ; // tentative shortest path length
100 3955 : GrB_Vector tmasked = NULL ;
101 3955 : GrB_Vector tReq = NULL ;
102 3955 : GrB_Vector tless = NULL ;
103 3955 : GrB_Vector s = NULL ;
104 3955 : GrB_Vector reach = NULL ;
105 3955 : GrB_Vector Empty = NULL ;
106 :
107 3955 : LG_TRY (LAGraph_CheckGraph (G, msg)) ;
108 3955 : LG_ASSERT (path_length != NULL && Delta != NULL, GrB_NULL_POINTER) ;
109 3955 : (*path_length) = NULL ;
110 :
111 : GrB_Index nvals ;
112 3955 : LG_TRY (GrB_Scalar_nvals (&nvals, Delta)) ;
113 3955 : LG_ASSERT_MSG (nvals == 1, GrB_EMPTY_OBJECT, "Delta is missing") ;
114 :
115 3954 : GrB_Matrix A = G->A ;
116 : GrB_Index n ;
117 3954 : GRB_TRY (GrB_Matrix_nrows (&n, A)) ;
118 3954 : 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 3954 : LG_TRY (LAGraph_Matrix_TypeName (typename, A, msg)) ;
128 3954 : LG_TRY (LAGraph_TypeFromName (&etype, typename, msg)) ;
129 :
130 3954 : GRB_TRY (GrB_Scalar_new (&lBound, etype)) ;
131 3926 : GRB_TRY (GrB_Scalar_new (&uBound, etype)) ;
132 3898 : GRB_TRY (GrB_Vector_new (&t, etype, n)) ;
133 3870 : GRB_TRY (GrB_Vector_new (&tmasked, etype, n)) ;
134 3842 : GRB_TRY (GrB_Vector_new (&tReq, etype, n)) ;
135 3814 : GRB_TRY (GrB_Vector_new (&Empty, GrB_BOOL, n)) ;
136 3786 : GRB_TRY (GrB_Vector_new (&tless, GrB_BOOL, n)) ;
137 3758 : GRB_TRY (GrB_Vector_new (&s, GrB_BOOL, n)) ;
138 3730 : GRB_TRY (GrB_Vector_new (&reach, GrB_BOOL, n)) ;
139 :
140 : #if LAGRAPH_SUITESPARSE
141 : // optional hints for SuiteSparse:GraphBLAS
142 3702 : GRB_TRY (GxB_set (t, GxB_SPARSITY_CONTROL, GxB_BITMAP)) ;
143 3674 : GRB_TRY (GxB_set (tmasked, GxB_SPARSITY_CONTROL, GxB_SPARSE)) ;
144 3674 : GRB_TRY (GxB_set (tReq, GxB_SPARSITY_CONTROL, GxB_SPARSE)) ;
145 3674 : GRB_TRY (GxB_set (tless, GxB_SPARSITY_CONTROL, GxB_SPARSE)) ;
146 3674 : GRB_TRY (GxB_set (s, GxB_SPARSITY_CONTROL, GxB_SPARSE)) ;
147 3674 : GRB_TRY (GxB_set (reach, GxB_SPARSITY_CONTROL, GxB_BITMAP)) ;
148 : #endif
149 :
150 : // select the operators, and set t (:) = infinity
151 : GrB_IndexUnaryOp ne, le, ge, lt, gt ;
152 : GrB_BinaryOp less_than ;
153 : GrB_Semiring min_plus ;
154 : int tcode ;
155 : int32_t delta_int32 ;
156 : int64_t delta_int64 ;
157 : uint32_t delta_uint32 ;
158 : uint64_t delta_uint64 ;
159 : float delta_fp32 ;
160 : double delta_fp64 ;
161 :
162 3646 : bool negative_edge_weights = true ;
163 :
164 3646 : if (etype == GrB_INT32)
165 : {
166 3419 : GRB_TRY (GrB_Scalar_extractElement (&delta_int32, Delta)) ;
167 3419 : GRB_TRY (GrB_assign (t, NULL, NULL, (int32_t) INT32_MAX,
168 : GrB_ALL, n, NULL)) ;
169 3391 : ne = GrB_VALUENE_INT32 ;
170 3391 : le = GrB_VALUELE_INT32 ;
171 3391 : ge = GrB_VALUEGE_INT32 ;
172 3391 : lt = GrB_VALUELT_INT32 ;
173 3391 : gt = GrB_VALUEGT_INT32 ;
174 3391 : less_than = GrB_LT_INT32 ;
175 3391 : min_plus = GrB_MIN_PLUS_SEMIRING_INT32 ;
176 3391 : tcode = 0 ;
177 : }
178 227 : else if (etype == GrB_INT64)
179 : {
180 48 : GRB_TRY (GrB_Scalar_extractElement (&delta_int64, Delta)) ;
181 48 : GRB_TRY (GrB_assign (t, NULL, NULL, (int64_t) INT64_MAX,
182 : GrB_ALL, n, NULL)) ;
183 48 : ne = GrB_VALUENE_INT64 ;
184 48 : le = GrB_VALUELE_INT64 ;
185 48 : ge = GrB_VALUEGE_INT64 ;
186 48 : lt = GrB_VALUELT_INT64 ;
187 48 : gt = GrB_VALUEGT_INT64 ;
188 48 : less_than = GrB_LT_INT64 ;
189 48 : min_plus = GrB_MIN_PLUS_SEMIRING_INT64 ;
190 48 : tcode = 1 ;
191 : }
192 179 : else if (etype == GrB_UINT32)
193 : {
194 12 : GRB_TRY (GrB_Scalar_extractElement (&delta_uint32, Delta)) ;
195 12 : GRB_TRY (GrB_assign (t, NULL, NULL, (uint32_t) UINT32_MAX,
196 : GrB_ALL, n, NULL)) ;
197 12 : ne = GrB_VALUENE_UINT32 ;
198 12 : le = GrB_VALUELE_UINT32 ;
199 12 : ge = GrB_VALUEGE_UINT32 ;
200 12 : lt = GrB_VALUELT_UINT32 ;
201 12 : gt = GrB_VALUEGT_UINT32 ;
202 12 : less_than = GrB_LT_UINT32 ;
203 12 : min_plus = GrB_MIN_PLUS_SEMIRING_UINT32 ;
204 12 : tcode = 2 ;
205 12 : negative_edge_weights = false ;
206 : }
207 167 : else if (etype == GrB_UINT64)
208 : {
209 12 : GRB_TRY (GrB_Scalar_extractElement (&delta_uint64, Delta)) ;
210 12 : GRB_TRY (GrB_assign (t, NULL, NULL, (uint64_t) UINT64_MAX,
211 : GrB_ALL, n, NULL)) ;
212 12 : ne = GrB_VALUENE_UINT64 ;
213 12 : le = GrB_VALUELE_UINT64 ;
214 12 : ge = GrB_VALUEGE_UINT64 ;
215 12 : lt = GrB_VALUELT_UINT64 ;
216 12 : gt = GrB_VALUEGT_UINT64 ;
217 12 : less_than = GrB_LT_UINT64 ;
218 12 : min_plus = GrB_MIN_PLUS_SEMIRING_UINT64 ;
219 12 : tcode = 3 ;
220 12 : negative_edge_weights = false ;
221 : }
222 155 : else if (etype == GrB_FP32)
223 : {
224 9 : GRB_TRY (GrB_Scalar_extractElement (&delta_fp32, Delta)) ;
225 9 : GRB_TRY (GrB_assign (t, NULL, NULL, (float) INFINITY,
226 : GrB_ALL, n, NULL)) ;
227 9 : ne = GrB_VALUENE_FP32 ;
228 9 : le = GrB_VALUELE_FP32 ;
229 9 : ge = GrB_VALUEGE_FP32 ;
230 9 : lt = GrB_VALUELT_FP32 ;
231 9 : gt = GrB_VALUEGT_FP32 ;
232 9 : less_than = GrB_LT_FP32 ;
233 9 : min_plus = GrB_MIN_PLUS_SEMIRING_FP32 ;
234 9 : tcode = 4 ;
235 : }
236 146 : else if (etype == GrB_FP64)
237 : {
238 145 : GRB_TRY (GrB_Scalar_extractElement (&delta_fp64, Delta)) ;
239 145 : GRB_TRY (GrB_assign (t, NULL, NULL, (double) INFINITY,
240 : GrB_ALL, n, NULL)) ;
241 145 : ne = GrB_VALUENE_FP64 ;
242 145 : le = GrB_VALUELE_FP64 ;
243 145 : ge = GrB_VALUEGE_FP64 ;
244 145 : lt = GrB_VALUELT_FP64 ;
245 145 : gt = GrB_VALUEGT_FP64 ;
246 145 : less_than = GrB_LT_FP64 ;
247 145 : min_plus = GrB_MIN_PLUS_SEMIRING_FP64 ;
248 145 : tcode = 5 ;
249 : }
250 : else
251 : {
252 1 : LG_ASSERT_MSG (false, GrB_NOT_IMPLEMENTED, "type not supported") ;
253 : }
254 :
255 : // check if the graph might have negative edge weights
256 3617 : if (negative_edge_weights)
257 : {
258 3593 : double emin = -1 ;
259 3593 : if (G->emin != NULL &&
260 452 : (G->emin_state == LAGraph_VALUE ||
261 238 : G->emin_state == LAGraph_BOUND))
262 : {
263 214 : GRB_TRY (GrB_Scalar_extractElement_FP64 (&emin, G->emin)) ;
264 : }
265 : // else
266 : // {
267 : // // a future Basic algorithm should compute G->emin and perhaps
268 : // // G->emax, then compute Delta automatically.
269 : // }
270 3593 : negative_edge_weights = (emin < 0) ;
271 : }
272 :
273 : // t (src) = 0
274 3617 : GRB_TRY (GrB_Vector_setElement (t, 0, source)) ;
275 :
276 : // reach (src) = true
277 3603 : GRB_TRY (GrB_Vector_setElement (reach, true, source)) ;
278 :
279 : // s (src) = true
280 3589 : GRB_TRY (GrB_Vector_setElement (s, true, source)) ;
281 :
282 : // AL = A .* (A <= Delta)
283 3547 : GRB_TRY (GrB_Matrix_new (&AL, etype, n, n)) ;
284 3505 : GRB_TRY (GrB_select (AL, NULL, NULL, le, A, Delta, NULL)) ;
285 3454 : GRB_TRY (GrB_wait (AL, GrB_MATERIALIZE)) ;
286 :
287 : // FUTURE: costly for some problems, taking up to 50% of the total time:
288 : // AH = A .* (A > Delta)
289 3454 : GRB_TRY (GrB_Matrix_new (&AH, etype, n, n)) ;
290 3412 : GRB_TRY (GrB_select (AH, NULL, NULL, gt, A, Delta, NULL)) ;
291 3336 : GRB_TRY (GrB_wait (AH, GrB_MATERIALIZE)) ;
292 :
293 : //--------------------------------------------------------------------------
294 : // while (t >= step*Delta) not empty
295 : //--------------------------------------------------------------------------
296 :
297 3336 : for (int64_t step = 0 ; ; step++)
298 9573 : {
299 :
300 : //----------------------------------------------------------------------
301 : // tmasked = all entries in t<reach> that are less than (step+1)*Delta
302 : //----------------------------------------------------------------------
303 :
304 12909 : setelement (uBound, (step+1)) ; // uBound = (step+1) * Delta
305 15596 : GRB_TRY (GrB_Vector_clear (tmasked)) ;
306 :
307 : // tmasked<reach> = t
308 : // FUTURE: this is costly, typically using Method 06s in SuiteSparse,
309 : // which is a very general-purpose one. Write a specialized kernel to
310 : // exploit the fact that reach and t are bitmap and tmasked starts
311 : // empty, or fuse this assignment with the GrB_select below.
312 12779 : GRB_TRY (GrB_assign (tmasked, reach, NULL, t, GrB_ALL, n, NULL)) ;
313 : // tmasked = select (tmasked < (step+1)*Delta)
314 12511 : GRB_TRY (GrB_select (tmasked, NULL, NULL, lt, tmasked, uBound, NULL)) ;
315 : // --- alternative:
316 : // FUTURE this is slower than the above but should be much faster.
317 : // GrB_select is computing a bitmap result then converting it to
318 : // sparse. t and reach are both bitmap and tmasked finally sparse.
319 : // tmasked<reach> = select (t < (step+1)*Delta)
320 : // GRB_TRY (GrB_select (tmasked, reach, NULL, lt, t, uBound, NULL)) ;
321 :
322 : GrB_Index tmasked_nvals ;
323 12343 : GRB_TRY (GrB_Vector_nvals (&tmasked_nvals, tmasked)) ;
324 :
325 : //----------------------------------------------------------------------
326 : // continue while the current bucket (tmasked) is not empty
327 : //----------------------------------------------------------------------
328 :
329 25660 : while (tmasked_nvals > 0)
330 : {
331 : // tReq = AL'*tmasked using the min_plus semiring
332 20033 : GRB_TRY (GrB_vxm (tReq, NULL, NULL, min_plus, tmasked, AL, NULL)) ;
333 :
334 : // s<struct(tmasked)> = true
335 18360 : GRB_TRY (GrB_assign (s, tmasked, NULL, (bool) true, GrB_ALL, n,
336 : GrB_DESC_S)) ;
337 :
338 : // if nvals (tReq) is 0, no need to continue the rest of this loop
339 : GrB_Index tReq_nvals ;
340 18153 : GRB_TRY (GrB_Vector_nvals (&tReq_nvals, tReq)) ;
341 20223 : if (tReq_nvals == 0) break ;
342 :
343 : // tless = (tReq .< t) using set intersection
344 16523 : GRB_TRY (GrB_eWiseMult (tless, NULL, NULL, less_than, tReq, t,
345 : NULL)) ;
346 :
347 : // remove explicit zeros from tless so it can be used as a
348 : // structural mask
349 : GrB_Index tless_nvals ;
350 16310 : GRB_TRY (GrB_select (tless, NULL, NULL, ne, tless, 0, NULL)) ;
351 16103 : GRB_TRY (GrB_Vector_nvals (&tless_nvals, tless)) ;
352 16103 : if (tless_nvals == 0) break ;
353 :
354 : // update reachable node list/mask
355 : // reach<struct(tless)> = true
356 14033 : GRB_TRY (GrB_assign (reach, tless, NULL, (bool) true, GrB_ALL, n,
357 : GrB_DESC_S)) ;
358 :
359 : // tmasked<struct(tless)> = select (tReq < (step+1)*Delta)
360 14033 : GRB_TRY (GrB_Vector_clear (tmasked)) ;
361 13985 : GRB_TRY (GrB_select (tmasked, tless, NULL, lt, tReq, uBound,
362 : GrB_DESC_S)) ;
363 :
364 : // For general graph with some negative weights:
365 13581 : if (negative_edge_weights)
366 : {
367 : // If all entries of the graph are known to be positive, and
368 : // the entries of tmasked are at least step*Delta, tReq =
369 : // tmasked min.+ AL must be >= step*Delta. Therefore, there is
370 : // no need to perform this GrB_select with ge to find tmasked
371 : // >= step*Delta from tReq.
372 10710 : setelement (lBound, (step)) ; // lBound = step*Delta
373 : // tmasked = select entries in tmasked that are >= step*Delta
374 10710 : GRB_TRY (GrB_select (tmasked, NULL, NULL, ge, tmasked, lBound,
375 : NULL)) ;
376 : }
377 :
378 : // t<struct(tless)> = tReq
379 13317 : GRB_TRY (GrB_assign (t, tless, NULL, tReq, GrB_ALL, n, GrB_DESC_S));
380 13317 : GRB_TRY (GrB_Vector_nvals (&tmasked_nvals, tmasked)) ;
381 : }
382 :
383 : // tmasked<s> = t
384 10670 : GRB_TRY (GrB_Vector_clear (tmasked)) ;
385 10638 : GRB_TRY (GrB_assign (tmasked, s, NULL, t, GrB_ALL, n, GrB_DESC_S)) ;
386 :
387 : // tReq = AH'*tmasked using the min_plus semiring
388 10404 : GRB_TRY (GrB_vxm (tReq, NULL, NULL, min_plus, tmasked, AH, NULL)) ;
389 :
390 : // tless = (tReq .< t) using set intersection
391 10212 : GRB_TRY (GrB_eWiseMult (tless, NULL, NULL, less_than, tReq, t, NULL)) ;
392 :
393 : // t<tless> = tReq, which computes t = min (t, tReq)
394 10116 : GRB_TRY (GrB_assign (t, tless, NULL, tReq, GrB_ALL, n, NULL)) ;
395 :
396 : //----------------------------------------------------------------------
397 : // find out how many left to be computed
398 : //----------------------------------------------------------------------
399 :
400 : // update reachable node list
401 : // reach<tless> = true
402 10116 : GRB_TRY (GrB_assign (reach, tless, NULL, (bool) true, GrB_ALL, n,
403 : NULL)) ;
404 :
405 : // remove previous buckets
406 : // reach<struct(s)> = Empty
407 10116 : GRB_TRY (GrB_assign (reach, s, NULL, Empty, GrB_ALL, n, GrB_DESC_S)) ;
408 : GrB_Index nreach ;
409 10110 : GRB_TRY (GrB_Vector_nvals (&nreach, reach)) ;
410 10110 : if (nreach == 0) break ;
411 :
412 9591 : GRB_TRY (GrB_Vector_clear (s)) ; // clear s for the next iteration
413 : }
414 :
415 : //--------------------------------------------------------------------------
416 : // free workspace and return result
417 : //--------------------------------------------------------------------------
418 :
419 519 : (*path_length) = t ;
420 519 : LG_FREE_WORK ;
421 519 : return (GrB_SUCCESS) ;
422 : }
|
