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