Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LAGraph_BF_full.c: Bellman-Ford single-source shortest paths, returns tree
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 and Timothy A. Davis, Texas A&M University
15 :
16 : //------------------------------------------------------------------------------
17 :
18 : // LAGraph_BF_full_mxv: Bellman-Ford single source shortest paths, returning
19 : // both the path lengths and the shortest-path tree.
20 :
21 : // LAGraph_BF_full_mxv performs a Bellman-Ford to find out shortest path
22 : // length, parent nodes along the path and the hops (number of edges) in the
23 : // path from given source vertex s in the range of [0, n) on graph with n nodes.
24 : // It works almost the same as LAGraph_BF_full except that it performs update
25 : // using GrB_mxv instead of GrB_vxm, therefore, it require the input matrix as
26 : // the transpose of adjacency matrix A with size n by n. That is, the input
27 : // sparse matrix has entry AT(i, j) if there is edge from vertex j to vertex i
28 : // with weight w, then AT(i, j) = w. While same as LAGraph_BF_full, it requires
29 : // AT(i, i) = 0 for all 0 <= i < n.
30 :
31 : // LAGraph_BF_full_mxv returns GrB_SUCCESS if it succeeds. In this case, there
32 : // are no negative-weight cycles in the graph, and d, pi, and h are returned.
33 : // The vector d has d(k) as the shortest distance from s to k. pi(k) = p+1,
34 : // where p is the parent node of k-th node in the shortest path. In particular,
35 : // pi(s) = 0. h(k) = hop(s, k), the number of edges from s to k in the shortest
36 : // path.
37 :
38 : // If the graph has a negative-weight cycle, GrB_NO_VALUE is returned, and the
39 : // GrB_Vectors d(k), pi(k) and h(k) (i.e., *pd_output, *ppi_output and
40 : // *ph_output respectively) will be NULL when negative-weight cycle detected.
41 :
42 : // Otherwise, other errors such as GrB_OUT_OF_MEMORY, GrB_INVALID_OBJECT, and
43 : // so on, can be returned, if these errors are found by the underlying
44 : // GrB_* functions.
45 :
46 :
47 : //------------------------------------------------------------------------------
48 :
49 : #define LG_FREE_WORK \
50 : { \
51 : GrB_free(&d); \
52 : GrB_free(&dtmp); \
53 : GrB_free(&Atmp); \
54 : GrB_free(&BF_Tuple3); \
55 : GrB_free(&BF_lMIN_Tuple3); \
56 : GrB_free(&BF_PLUSrhs_Tuple3); \
57 : GrB_free(&BF_EQ_Tuple3); \
58 : GrB_free(&BF_lMIN_Tuple3_Monoid); \
59 : GrB_free(&BF_lMIN_PLUSrhs_Tuple3); \
60 : LAGraph_Free ((void**)&I, NULL); \
61 : LAGraph_Free ((void**)&J, NULL); \
62 : LAGraph_Free ((void**)&w, NULL); \
63 : LAGraph_Free ((void**)&W, NULL); \
64 : LAGraph_Free ((void**)&h, NULL); \
65 : LAGraph_Free ((void**)&pi, NULL); \
66 : }
67 :
68 : #define LG_FREE_ALL \
69 : { \
70 : LG_FREE_WORK ; \
71 : GrB_free (pd_output); \
72 : GrB_free (ppi_output); \
73 : GrB_free (ph_output); \
74 : }
75 :
76 : #include <LAGraph.h>
77 : #include <LAGraphX.h>
78 : #include <LG_internal.h> // from src/utility
79 :
80 : typedef void (*LAGraph_binary_function) (void *, const void *, const void *) ;
81 :
82 : //------------------------------------------------------------------------------
83 : // data type for each entry of the adjacent matrix A and "distance" vector d;
84 : // <INFINITY,INFINITY,INFINITY> corresponds to nonexistence of a path, and
85 : // the value <0, 0, NULL> corresponds to a path from a vertex to itself
86 : //------------------------------------------------------------------------------
87 :
88 : typedef struct
89 : {
90 : double w; // w corresponds to a path weight.
91 : GrB_Index h; // h corresponds to a path size or number of hops.
92 : GrB_Index pi;// pi corresponds to the penultimate vertex along a path.
93 : // vertex indexed as 1, 2, 3, ... , V, and pi = 0 (as nil)
94 : // for u=v, and pi = UINT64_MAX (as inf) for (u,v) not in E
95 : }
96 : BF_Tuple3_struct;
97 :
98 : //------------------------------------------------------------------------------
99 : // binary functions, z=f(x,y), where Tuple3xTuple3 -> Tuple3
100 : //------------------------------------------------------------------------------
101 :
102 20297 : void BF_lMIN_mxv
103 : (
104 : BF_Tuple3_struct *z,
105 : const BF_Tuple3_struct *y,
106 : const BF_Tuple3_struct *x
107 : )
108 : {
109 20297 : if (x->w < y->w
110 15239 : || (x->w == y->w && x->h < y->h)
111 15239 : || (x->w == y->w && x->h == y->h && x->pi < y->pi))
112 : {
113 5058 : if (z != x) { *z = *x; }
114 : }
115 : else
116 : {
117 15239 : *z = *y;
118 : }
119 20297 : }
120 :
121 25194 : void BF_PLUSrhs_mxv
122 : (
123 : BF_Tuple3_struct *z,
124 : const BF_Tuple3_struct *y,
125 : const BF_Tuple3_struct *x
126 : )
127 : {
128 25194 : z->w = x->w + y->w;
129 25194 : z->h = x->h + y->h;
130 25194 : z->pi = (x->pi != UINT64_MAX && y->pi != 0) ? y->pi : x->pi ;
131 25194 : }
132 :
133 4424 : void BF_EQ_mxv
134 : (
135 : bool *z,
136 : const BF_Tuple3_struct *y,
137 : const BF_Tuple3_struct *x
138 : )
139 : {
140 4424 : (*z) = (x->w == y->w && x->h == y->h && x->pi == y->pi) ;
141 4424 : }
142 :
143 : // Given the transpose of a n-by-n adjacency matrix A and a source vertex s.
144 : // If there is no negative-weight cycle reachable from s, return the distances
145 : // of shortest paths from s and parents along the paths as vector d. Otherwise,
146 : // returns d=NULL if there is a negtive-weight cycle.
147 : // pd_output is pointer to a GrB_Vector, where the i-th entry is d(s,i), the
148 : // sum of edges length in the shortest path
149 : // ppi_output is pointer to a GrB_Vector, where the i-th entry is pi(i), the
150 : // parent of i-th vertex in the shortest path
151 : // ph_output is pointer to a GrB_Vector, where the i-th entry is h(s,i), the
152 : // number of edges from s to i in the shortest path
153 : // AT has zeros on diagonal and weights on corresponding entries of edges
154 : // s is given index for source vertex
155 :
156 5 : GrB_Info LAGraph_BF_full_mxv
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 AT, //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, dtmp = 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_EQ_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 AT
181 5 : GrB_Index *h = NULL, *pi = NULL;
182 5 : double *w = NULL;
183 5 : BF_Tuple3_struct *W = NULL;
184 :
185 5 : LG_ASSERT (AT != NULL && pd_output != NULL &&
186 : ppi_output != NULL && ph_output != NULL, GrB_NULL_POINTER) ;
187 :
188 5 : *pd_output = NULL;
189 5 : *ppi_output = NULL;
190 5 : *ph_output = NULL;
191 5 : GRB_TRY (GrB_Matrix_nrows (&nrows, AT)) ;
192 5 : GRB_TRY (GrB_Matrix_ncols (&ncols, AT)) ;
193 5 : GRB_TRY (GrB_Matrix_nvals (&nz, AT));
194 5 : LG_ASSERT_MSG (nrows == ncols, -1002, "A must be square") ;
195 5 : n = nrows;
196 5 : LG_ASSERT_MSG (s < n, GrB_INVALID_INDEX, "invalid source node") ;
197 :
198 : //--------------------------------------------------------------------------
199 : // create all GrB_Type GrB_BinaryOp GrB_Monoid and GrB_Semiring
200 : //--------------------------------------------------------------------------
201 : // GrB_Type
202 5 : GRB_TRY (GrB_Type_new(&BF_Tuple3, sizeof(BF_Tuple3_struct)));
203 :
204 : // GrB_BinaryOp
205 5 : GRB_TRY (GrB_BinaryOp_new(&BF_EQ_Tuple3,
206 : (LAGraph_binary_function) (&BF_EQ_mxv),
207 : GrB_BOOL, BF_Tuple3, BF_Tuple3));
208 5 : GRB_TRY (GrB_BinaryOp_new(&BF_lMIN_Tuple3,
209 : (LAGraph_binary_function) (&BF_lMIN_mxv),
210 : BF_Tuple3, BF_Tuple3, BF_Tuple3));
211 5 : GRB_TRY (GrB_BinaryOp_new(&BF_PLUSrhs_Tuple3,
212 : (LAGraph_binary_function)(&BF_PLUSrhs_mxv),
213 : BF_Tuple3, BF_Tuple3, BF_Tuple3));
214 :
215 : // GrB_Monoid
216 5 : BF_Tuple3_struct BF_identity = (BF_Tuple3_struct) { .w = INFINITY,
217 : .h = UINT64_MAX, .pi = UINT64_MAX };
218 5 : GRB_TRY (GrB_Monoid_new_UDT(&BF_lMIN_Tuple3_Monoid, BF_lMIN_Tuple3,
219 : &BF_identity));
220 :
221 : //GrB_Semiring
222 5 : GRB_TRY (GrB_Semiring_new(&BF_lMIN_PLUSrhs_Tuple3,
223 : BF_lMIN_Tuple3_Monoid, BF_PLUSrhs_Tuple3));
224 :
225 : //--------------------------------------------------------------------------
226 : // allocate arrays used for tuplets
227 : //--------------------------------------------------------------------------
228 :
229 5 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &I, nz, sizeof(GrB_Index), msg)) ;
230 5 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &J, nz, sizeof(GrB_Index), msg)) ;
231 5 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &w, nz, sizeof(double), msg)) ;
232 5 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &W, nz, sizeof(BF_Tuple3_struct),
233 : msg)) ;
234 :
235 : //--------------------------------------------------------------------------
236 : // create matrix Atmp based on AT, while its entries become BF_Tuple3 type
237 : //--------------------------------------------------------------------------
238 5 : GRB_TRY (GrB_Matrix_extractTuples_FP64(I, J, w, &nz, AT));
239 951 : for (GrB_Index k = 0; k < nz; k++)
240 : {
241 946 : if (w[k] == 0) //diagonal entries
242 : {
243 184 : W[k] = (BF_Tuple3_struct) { .w = 0, .h = 0, .pi = 0 };
244 : }
245 : else
246 : {
247 : // w(i, j) is the weight of edge from j to i, so the parent here
248 : // should be j
249 762 : W[k] = (BF_Tuple3_struct) { .w = w[k], .h = 1, .pi = J[k] + 1 };
250 : }
251 : }
252 5 : GRB_TRY (GrB_Matrix_new(&Atmp, BF_Tuple3, n, n));
253 5 : GRB_TRY (GrB_Matrix_build_UDT(Atmp, I, J, W, nz, BF_lMIN_Tuple3));
254 5 : LAGraph_Free ((void**)&I, NULL);
255 5 : LAGraph_Free ((void**)&J, NULL);
256 5 : LAGraph_Free ((void**)&W, NULL);
257 5 : LAGraph_Free ((void**)&w, NULL);
258 :
259 : //--------------------------------------------------------------------------
260 : // create and initialize "distance" vector d
261 : //--------------------------------------------------------------------------
262 5 : GRB_TRY (GrB_Vector_new(&d, BF_Tuple3, n));
263 : // initial distance from s to itself
264 5 : BF_Tuple3_struct d0 = (BF_Tuple3_struct) { .w = 0, .h = 0, .pi = 0 };
265 5 : GRB_TRY (GrB_Vector_setElement_UDT(d, &d0, s));
266 :
267 : //--------------------------------------------------------------------------
268 : // start the Bellman Ford process
269 : //--------------------------------------------------------------------------
270 : // copy d to dtmp in order to create a same size of vector
271 5 : GRB_TRY (GrB_Vector_dup(&dtmp, d));
272 5 : bool same= false; // variable indicating if d == dtmp
273 5 : int64_t iter = 0; // number of iterations
274 :
275 : // terminate when no new path is found or more than V-1 loops
276 93 : while (!same && iter < n - 1)
277 : {
278 : // execute semiring on d and AT, and save the result to dtmp
279 88 : GRB_TRY (GrB_mxv(dtmp, GrB_NULL, GrB_NULL, BF_lMIN_PLUSrhs_Tuple3,
280 : Atmp, d, GrB_NULL));
281 88 : LG_TRY (LAGraph_Vector_IsEqualOp (&same, dtmp, d, BF_EQ_Tuple3, NULL));
282 88 : if (!same)
283 : {
284 85 : GrB_Vector ttmp = dtmp;
285 85 : dtmp = d;
286 85 : d = ttmp;
287 : }
288 88 : iter ++;
289 : }
290 :
291 : // check for negative-weight cycle only when there was a new path in the
292 : // last loop, otherwise, there can't be a negative-weight cycle.
293 5 : if (!same)
294 : {
295 : // execute semiring again to check for negative-weight cycle
296 2 : GRB_TRY (GrB_mxv(dtmp, GrB_NULL, GrB_NULL, BF_lMIN_PLUSrhs_Tuple3,
297 : Atmp, d, GrB_NULL));
298 2 : LG_TRY (LAGraph_Vector_IsEqualOp (&same, dtmp, d, BF_EQ_Tuple3, NULL));
299 :
300 : // if d != dtmp, then there is a negative-weight cycle in the graph
301 2 : if (!same)
302 : {
303 : // printf("A negative-weight cycle found. \n");
304 2 : LG_FREE_ALL;
305 2 : return (GrB_NO_VALUE) ;
306 : }
307 : }
308 :
309 : //--------------------------------------------------------------------------
310 : // extract tuple from "distance" vector d and create GrB_Vectors for output
311 : //--------------------------------------------------------------------------
312 :
313 3 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &I, n, sizeof(GrB_Index), msg)) ;
314 3 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &W, n, sizeof(BF_Tuple3_struct),
315 : msg)) ;
316 3 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &w, n, sizeof(double), msg)) ;
317 3 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &h, n, sizeof(GrB_Index), msg)) ;
318 3 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &pi, n, sizeof(GrB_Index), msg)) ;
319 :
320 3 : nz = n ;
321 3 : GRB_TRY (GrB_Vector_extractTuples_UDT (I, (void *) W, &nz, d));
322 :
323 111 : for (GrB_Index k = 0; k < n; k++)
324 : {
325 108 : w [k] = W[k].w ;
326 108 : h [k] = W[k].h ;
327 108 : pi[k] = W[k].pi;
328 : }
329 :
330 3 : GRB_TRY (GrB_Vector_new(pd_output, GrB_FP64, n));
331 3 : GRB_TRY (GrB_Vector_new(ppi_output, GrB_UINT64, n));
332 3 : GRB_TRY (GrB_Vector_new(ph_output, GrB_UINT64, n));
333 3 : GRB_TRY (GrB_Vector_build (*pd_output , I, w , nz, GrB_MIN_FP64 ));
334 3 : GRB_TRY (GrB_Vector_build (*ppi_output, I, pi, nz, GrB_MIN_UINT64));
335 3 : GRB_TRY (GrB_Vector_build (*ph_output , I, h , nz, GrB_MIN_UINT64));
336 3 : LG_FREE_WORK;
337 3 : return (GrB_SUCCESS) ;
338 : }
|
