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