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 : }
|
