Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LAGraph_BF_full2.c: Bellman-Ford single-source shortest paths, returns tree,
3 : // while diagonal of input matrix A needs not to be explicit 0, using the
4 : // frontier idea from Roi Lipman
5 : //------------------------------------------------------------------------------
6 :
7 : // LAGraph, (c) 2019-2022 by The LAGraph Contributors, All Rights Reserved.
8 : // SPDX-License-Identifier: BSD-2-Clause
9 : //
10 : // For additional details (including references to third party source code and
11 : // other files) see the LICENSE file or contact permission@sei.cmu.edu. See
12 : // Contributors.txt for a full list of contributors. Created, in part, with
13 : // funding and support from the U.S. Government (see Acknowledgments.txt file).
14 : // DM22-0790
15 :
16 : // Contributed by Jinhao Chen and Timothy A. Davis, Texas A&M University
17 :
18 : //------------------------------------------------------------------------------
19 :
20 : // LAGraph_BF_full2: Bellman-Ford single source shortest paths, returning both
21 : // the path lengths and the shortest-path tree.
22 :
23 : // LAGraph_BF_full2 performs a Bellman-Ford to find out shortest path, parent
24 : // nodes along the path and the hops (number of edges) in the path from given
25 : // source vertex s in the range of [0, n) on graph given as matrix A with size
26 : // n*n. The sparse matrix A has entry A(i, j) if there is an edge from vertex i
27 : // to vertex j with weight w, then A(i, j) = w.
28 :
29 : // LAGraph_BF_full2 returns GrB_SUCCESS if it succeeds. In this case, there
30 : // are no negative-weight cycles in the graph, and d, pi, and h are returned.
31 : // The vector d has d(k) as the shortest distance from s to k. pi(k) = p+1,
32 : // where p is the parent node of k-th node in the shortest path. In particular,
33 : // pi(s) = 0. h(k) = hop(s, k), the number of edges from s to k in the shortest
34 : // path.
35 :
36 : // If the graph has a negative-weight cycle, GrB_NO_VALUE is returned, and the
37 : // GrB_Vectors d(k), pi(k) and h(k) (i.e., *pd_output, *ppi_output and
38 : // *ph_output respectively) will be NULL when negative-weight cycle detected.
39 :
40 : // Otherwise, other errors such as GrB_OUT_OF_MEMORY, GrB_INVALID_OBJECT, and
41 : // so on, can be returned, if these errors are found by the underlying
42 : // GrB_* functions.
43 :
44 : //------------------------------------------------------------------------------
45 :
46 : #define LG_FREE_WORK \
47 : { \
48 : GrB_free(&d); \
49 : GrB_free(&dtmp); \
50 : GrB_free(&dfrontier); \
51 : GrB_free(&Atmp); \
52 : GrB_free(&BF_Tuple3); \
53 : GrB_free(&BF_lMIN_Tuple3); \
54 : GrB_free(&BF_PLUSrhs_Tuple3); \
55 : GrB_free(&BF_EQ_Tuple3); \
56 : GrB_free(&BF_lMIN_Tuple3_Monoid); \
57 : GrB_free(&BF_lMIN_PLUSrhs_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 : BF2_Tuple3_struct;
95 :
96 : //------------------------------------------------------------------------------
97 : // binary functions, z=f(x,y), where Tuple3xTuple3 -> Tuple3
98 : //------------------------------------------------------------------------------
99 :
100 20299 : void BF2_lMIN2
101 : (
102 : BF2_Tuple3_struct *z,
103 : const BF2_Tuple3_struct *x,
104 : const BF2_Tuple3_struct *y
105 : )
106 : {
107 20299 : if (x->w < y->w
108 8956 : || (x->w == y->w && x->h < y->h)
109 8956 : || (x->w == y->w && x->h == y->h && x->pi < y->pi))
110 : {
111 11626 : if (z != x) { *z = *x; }
112 : }
113 : else
114 : {
115 8673 : *z = *y;
116 : }
117 20299 : }
118 :
119 20476 : void BF2_PLUSrhs2
120 : (
121 : BF2_Tuple3_struct *z,
122 : const BF2_Tuple3_struct *x,
123 : const BF2_Tuple3_struct *y
124 : )
125 : {
126 20476 : z->w = x->w + y->w ;
127 20476 : z->h = x->h + y->h ;
128 20476 : z->pi = (x->pi != UINT64_MAX && y->pi != 0) ? y->pi : x->pi ;
129 20476 : }
130 :
131 4424 : void BF2_EQ
132 : (
133 : bool *z,
134 : const BF2_Tuple3_struct *x,
135 : const BF2_Tuple3_struct *y
136 : )
137 : {
138 4424 : (*z) = (x->w == y->w && x->h == y->h && x->pi == y->pi) ;
139 4424 : }
140 :
141 : // Given a n-by-n adjacency matrix A and a source vertex s.
142 : // If there is no negative-weight cycle reachable from s, return the distances
143 : // of shortest paths from s and parents along the paths as vector d. Otherwise,
144 : // returns d=NULL if there is a negtive-weight cycle.
145 : // pd_output is pointer to a GrB_Vector, where the i-th entry is d(s,i), the
146 : // sum of edges length in the shortest path
147 : // ppi_output is pointer to a GrB_Vector, where the i-th entry is pi(i), the
148 : // parent of i-th vertex in the shortest path
149 : // ph_output is pointer to a GrB_Vector, where the i-th entry is h(s,i), the
150 : // number of edges from s to i in the shortest path
151 : // A has weights on corresponding entries of edges
152 : // s is given index for source vertex
153 :
154 5 : GrB_Info LAGraph_BF_full2
155 : (
156 : GrB_Vector *pd_output, //the pointer to the vector of distance
157 : GrB_Vector *ppi_output, //the pointer to the vector of parent
158 : GrB_Vector *ph_output, //the pointer to the vector of hops
159 : const GrB_Matrix A, //matrix for the graph
160 : const GrB_Index s //given index of the source
161 : )
162 : {
163 : GrB_Info info;
164 5 : char *msg = NULL ;
165 : // tmp vector to store distance vector after n (i.e., V) loops
166 5 : GrB_Vector d = NULL, dtmp = NULL, dfrontier = NULL;
167 5 : GrB_Matrix Atmp = NULL;
168 : GrB_Type BF_Tuple3;
169 :
170 : GrB_BinaryOp BF_lMIN_Tuple3;
171 : GrB_BinaryOp BF_PLUSrhs_Tuple3;
172 : GrB_BinaryOp BF_EQ_Tuple3;
173 :
174 : GrB_Monoid BF_lMIN_Tuple3_Monoid;
175 : GrB_Semiring BF_lMIN_PLUSrhs_Tuple3;
176 :
177 : GrB_Index nrows, ncols, n, nz; // n = # of row/col, nz = # of nnz in graph
178 5 : GrB_Index *I = NULL, *J = NULL; // for col/row indices of entries from A
179 5 : GrB_Index *h = NULL, *pi = NULL;
180 5 : double *w = NULL;
181 5 : BF2_Tuple3_struct *W = NULL;
182 :
183 5 : LG_ASSERT (A != NULL && pd_output != NULL &&
184 : ppi_output != NULL && ph_output != NULL, GrB_NULL_POINTER) ;
185 :
186 5 : *pd_output = NULL;
187 5 : *ppi_output = NULL;
188 5 : *ph_output = NULL;
189 5 : GRB_TRY (GrB_Matrix_nrows (&nrows, A)) ;
190 5 : GRB_TRY (GrB_Matrix_ncols (&ncols, A)) ;
191 5 : GRB_TRY (GrB_Matrix_nvals (&nz, A));
192 5 : LG_ASSERT_MSG (nrows == ncols, -1002, "A must be square") ;
193 5 : n = nrows;
194 5 : LG_ASSERT_MSG (s < n, GrB_INVALID_INDEX, "invalid source node") ;
195 :
196 : //--------------------------------------------------------------------------
197 : // create all GrB_Type GrB_BinaryOp GrB_Monoid and GrB_Semiring
198 : //--------------------------------------------------------------------------
199 : // GrB_Type
200 5 : GRB_TRY (GrB_Type_new(&BF_Tuple3, sizeof(BF2_Tuple3_struct)));
201 :
202 : // GrB_BinaryOp
203 5 : GRB_TRY (GrB_BinaryOp_new(&BF_EQ_Tuple3,
204 : (LAGraph_binary_function) (&BF2_EQ), GrB_BOOL, BF_Tuple3, BF_Tuple3));
205 5 : GRB_TRY (GrB_BinaryOp_new(&BF_lMIN_Tuple3,
206 : (LAGraph_binary_function) (&BF2_lMIN2),
207 : BF_Tuple3, BF_Tuple3, BF_Tuple3));
208 5 : GRB_TRY (GrB_BinaryOp_new(&BF_PLUSrhs_Tuple3,
209 : (LAGraph_binary_function)(&BF2_PLUSrhs2),
210 : BF_Tuple3, BF_Tuple3, BF_Tuple3));
211 :
212 : // GrB_Monoid
213 5 : BF2_Tuple3_struct BF_identity = (BF2_Tuple3_struct) { .w = INFINITY,
214 : .h = UINT64_MAX, .pi = UINT64_MAX };
215 5 : GRB_TRY (GrB_Monoid_new_UDT(&BF_lMIN_Tuple3_Monoid, BF_lMIN_Tuple3,
216 : &BF_identity));
217 :
218 : //GrB_Semiring
219 5 : GRB_TRY (GrB_Semiring_new(&BF_lMIN_PLUSrhs_Tuple3,
220 : BF_lMIN_Tuple3_Monoid, BF_PLUSrhs_Tuple3));
221 :
222 : //--------------------------------------------------------------------------
223 : // allocate arrays used for tuplets
224 : //--------------------------------------------------------------------------
225 :
226 5 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &I, nz, sizeof(GrB_Index), msg)) ;
227 5 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &J, nz, sizeof(GrB_Index), msg)) ;
228 5 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &w, nz, sizeof(double), msg)) ;
229 5 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &W, nz, sizeof(BF2_Tuple3_struct),
230 : msg)) ;
231 :
232 : //--------------------------------------------------------------------------
233 : // create matrix Atmp based on A, while its entries become BF_Tuple3 type
234 : //--------------------------------------------------------------------------
235 :
236 5 : GRB_TRY (GrB_Matrix_extractTuples_FP64(I, J, w, &nz, A));
237 : int nthreads, nthreads_outer, nthreads_inner ;
238 5 : LG_TRY (LAGraph_GetNumThreads (&nthreads_outer, &nthreads_inner, msg)) ;
239 5 : nthreads = nthreads_outer * nthreads_inner ;
240 5 : printf ("nthreads %d\n", nthreads) ;
241 : int64_t k;
242 : #pragma omp parallel for num_threads(nthreads) schedule(static)
243 773 : for (k = 0; k < nz; k++)
244 : {
245 768 : if (w[k] == 0) //diagonal entries
246 : {
247 2 : W[k] = (BF2_Tuple3_struct) { .w = 0, .h = 0, .pi = 0 };
248 : }
249 : else
250 : {
251 766 : W[k] = (BF2_Tuple3_struct) { .w = w[k], .h = 1, .pi = I[k] + 1 };
252 : }
253 : }
254 5 : GRB_TRY (GrB_Matrix_new(&Atmp, BF_Tuple3, n, n));
255 5 : GRB_TRY (GrB_Matrix_build_UDT(Atmp, I, J, W, nz, BF_lMIN_Tuple3));
256 5 : LAGraph_Free ((void**)&I, NULL);
257 5 : LAGraph_Free ((void**)&J, NULL);
258 5 : LAGraph_Free ((void**)&W, NULL);
259 5 : LAGraph_Free ((void**)&w, NULL);
260 :
261 : //--------------------------------------------------------------------------
262 : // create and initialize "distance" vector d
263 : //--------------------------------------------------------------------------
264 5 : GRB_TRY (GrB_Vector_new(&d, BF_Tuple3, n));
265 : // initial distance from s to itself
266 5 : BF2_Tuple3_struct d0 = (BF2_Tuple3_struct) { .w = 0, .h = 0, .pi = 0 };
267 5 : GRB_TRY (GrB_Vector_setElement_UDT(d, &d0, s));
268 :
269 : //--------------------------------------------------------------------------
270 : // start the Bellman Ford process
271 : //--------------------------------------------------------------------------
272 : // copy d to dtmp in order to create a same size of vector
273 5 : GRB_TRY (GrB_Vector_dup(&dtmp, d));
274 5 : GRB_TRY (GrB_Vector_dup(&dfrontier, d));
275 5 : bool same= false; // variable indicating if d == dtmp
276 5 : int64_t iter = 0; // number of iterations
277 :
278 : // terminate when no new path is found or more than V-1 loops
279 93 : while (!same && iter < n - 1)
280 : {
281 : // execute semiring on d and A, and save the result to dtmp
282 88 : GRB_TRY (GrB_vxm(dfrontier, GrB_NULL, GrB_NULL,
283 : BF_lMIN_PLUSrhs_Tuple3, dfrontier, Atmp, GrB_NULL));
284 :
285 : // dtmp[i] = min(d[i], dfrontier[i]).
286 88 : GrB_Vector_eWiseAdd_BinaryOp(dtmp, GrB_NULL, GrB_NULL, BF_lMIN_Tuple3,
287 : d, dfrontier, GrB_NULL);
288 :
289 88 : LG_TRY (LAGraph_Vector_IsEqualOp (&same, dtmp, d, BF_EQ_Tuple3, NULL));
290 88 : if (!same)
291 : {
292 85 : GrB_Vector ttmp = dtmp;
293 85 : dtmp = d;
294 85 : d = ttmp;
295 : }
296 88 : iter ++;
297 : }
298 :
299 : // check for negative-weight cycle only when there was a new path in the
300 : // last loop, otherwise, there can't be a negative-weight cycle.
301 5 : if (!same)
302 : {
303 : // execute semiring again to check for negative-weight cycle
304 2 : GRB_TRY (GrB_vxm(dfrontier, GrB_NULL, GrB_NULL,
305 : BF_lMIN_PLUSrhs_Tuple3, dfrontier, Atmp, GrB_NULL));
306 :
307 : // dtmp[i] = min(d[i], dfrontier[i]).
308 2 : GrB_Vector_eWiseAdd_BinaryOp(dtmp, GrB_NULL, GrB_NULL, BF_lMIN_Tuple3,
309 : d, dfrontier, GrB_NULL);
310 :
311 : // if d != dtmp, then there is a negative-weight cycle in the graph
312 2 : LG_TRY (LAGraph_Vector_IsEqualOp (&same, dtmp, d, BF_EQ_Tuple3, NULL));
313 2 : if (!same)
314 : {
315 : // printf("A negative-weight cycle found. \n");
316 2 : LG_FREE_ALL;
317 2 : return (GrB_NO_VALUE) ;
318 : }
319 : }
320 :
321 : //--------------------------------------------------------------------------
322 : // extract tuple from "distance" vector d and create GrB_Vectors for output
323 : //--------------------------------------------------------------------------
324 :
325 3 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &I, n, sizeof(GrB_Index), msg)) ;
326 3 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &W, n, sizeof(BF2_Tuple3_struct),
327 : msg)) ;
328 3 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &w, n, sizeof(double), msg)) ;
329 3 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &h, n, sizeof(GrB_Index), msg)) ;
330 3 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &pi, n, sizeof(GrB_Index), msg)) ;
331 :
332 3 : nz = n ;
333 3 : GRB_TRY (GrB_Vector_extractTuples_UDT (I, (void *) W, &nz, d));
334 :
335 111 : for (k = 0; k < nz; k++)
336 : {
337 108 : w [k] = W[k].w ;
338 108 : h [k] = W[k].h ;
339 108 : pi[k] = W[k].pi;
340 : }
341 3 : GRB_TRY (GrB_Vector_new(pd_output, GrB_FP64, n));
342 3 : GRB_TRY (GrB_Vector_new(ppi_output, GrB_UINT64, n));
343 3 : GRB_TRY (GrB_Vector_new(ph_output, GrB_UINT64, n));
344 3 : GRB_TRY (GrB_Vector_build (*pd_output , I, w , nz, GrB_MIN_FP64 ));
345 3 : GRB_TRY (GrB_Vector_build (*ppi_output, I, pi, nz, GrB_MIN_UINT64));
346 3 : GRB_TRY (GrB_Vector_build (*ph_output , I, h , nz, GrB_MIN_UINT64));
347 3 : LG_FREE_WORK;
348 3 : return (GrB_SUCCESS) ;
349 : }
|
