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: Bellman-Ford single source shortest paths, returning both
19 : // the path lengths and the shortest-path tree.
20 :
21 : // LAGraph_BF_full performs a Bellman-Ford to find out shortest path, parent
22 : // nodes along the path and the hops (number of edges) in the path from given
23 : // source vertex s in the range of [0, n) on graph given as matrix A with size
24 : // n*n. The sparse matrix A has entry A(i, j) if there is an edge from vertex i
25 : // to vertex j with weight w, then A(i, j) = w. Furthermore, LAGraph_BF_full
26 : // requires A(i, i) = 0 for all 0 <= i < n.
27 :
28 : // LAGraph_BF_full returns GrB_SUCCESS if successful, and GrB_NO_VALUE if it
29 : // detects the existence of negative- weight cycle. The GrB_Vector d(k), pi(k)
30 : // and h(k) (i.e., *pd_output, *ppi_output and *ph_output respectively) will
31 : // be NULL when negative-weight cycle detected. Otherwise, the vector d has
32 : // d(k) as the shortest distance from s to k. pi(k) = p+1, where p is the
33 : // parent node of k-th node in the shortest path. In particular, pi(s) = 0.
34 : // h(k) = hop(s, k), the number of edges from s to k in the shortest path.
35 :
36 : //------------------------------------------------------------------------------
37 :
38 : #define LG_FREE_ALL \
39 : { \
40 : GrB_free(&d); \
41 : GrB_free(&dtmp); \
42 : GrB_free(&Atmp); \
43 : GrB_free(&BF_Tuple3); \
44 : GrB_free(&BF_lMIN_Tuple3); \
45 : GrB_free(&BF_PLUSrhs_Tuple3); \
46 : GrB_free(&BF_EQ_Tuple3); \
47 : GrB_free(&BF_lMIN_Tuple3_Monoid); \
48 : GrB_free(&BF_lMIN_PLUSrhs_Tuple3); \
49 : LAGraph_Free ((void**)&I, NULL); \
50 : LAGraph_Free ((void**)&J, NULL); \
51 : LAGraph_Free ((void**)&w, NULL); \
52 : LAGraph_Free ((void**)&W, NULL); \
53 : LAGraph_Free ((void**)&h, NULL); \
54 : LAGraph_Free ((void**)&pi, NULL); \
55 : }
56 :
57 : #include <LAGraph.h>
58 : #include <LAGraphX.h>
59 : #include <LG_internal.h> // from src/utility
60 :
61 : typedef void (*LAGraph_binary_function) (void *, const void *, const void *) ;
62 :
63 : //------------------------------------------------------------------------------
64 : // data type for each entry of the adjacent matrix A and "distance" vector d;
65 : // <INFINITY,INFINITY,INFINITY> corresponds to nonexistence of a path, and
66 : // the value <0, 0, NULL> corresponds to a path from a vertex to itself
67 : //------------------------------------------------------------------------------
68 : typedef struct
69 : {
70 : double w; // w corresponds to a path weight.
71 : GrB_Index h; // h corresponds to a path size or number of hops.
72 : GrB_Index pi;// pi corresponds to the penultimate vertex along a path.
73 : // vertex indexed as 1, 2, 3, ... , V, and pi = 0 (as nil)
74 : // for u=v, and pi = UINT64_MAX (as inf) for (u,v) not in E
75 : }
76 : BF_Tuple3_struct;
77 :
78 : //------------------------------------------------------------------------------
79 : // 2 binary functions, z=f(x,y), where Tuple3xTuple3 -> Tuple3
80 : //------------------------------------------------------------------------------
81 :
82 20297 : void BF_lMIN
83 : (
84 : BF_Tuple3_struct *z,
85 : const BF_Tuple3_struct *x,
86 : const BF_Tuple3_struct *y
87 : )
88 : {
89 20297 : if (x->w < y->w
90 5788 : || (x->w == y->w && x->h < y->h)
91 5788 : || (x->w == y->w && x->h == y->h && x->pi < y->pi))
92 : {
93 14565 : if (z != x) { *z = *x; }
94 : }
95 : else
96 : {
97 5732 : *z = *y;
98 : }
99 20297 : }
100 :
101 25194 : void BF_PLUSrhs
102 : (
103 : BF_Tuple3_struct *z,
104 : const BF_Tuple3_struct *x,
105 : const BF_Tuple3_struct *y
106 : )
107 : {
108 25194 : z->w = x->w + y->w ;
109 25194 : z->h = x->h + y->h ;
110 25194 : z->pi = (x->pi != UINT64_MAX && y->pi != 0) ? y->pi : x->pi ;
111 25194 : }
112 :
113 4424 : void BF_EQ
114 : (
115 : bool *z,
116 : const BF_Tuple3_struct *x,
117 : const BF_Tuple3_struct *y
118 : )
119 : {
120 4424 : (*z) = (x->w == y->w && x->h == y->h && x->pi == y->pi) ;
121 4424 : }
122 :
123 : // Given a n-by-n adjacency matrix A and a source vertex s.
124 : // If there is no negative-weight cycle reachable from s, return the distances
125 : // of shortest paths from s and parents along the paths as vector d. Otherwise,
126 : // returns d=NULL if there is a negtive-weight cycle.
127 : // pd_output is pointer to a GrB_Vector, where the i-th entry is d(s,i), the
128 : // sum of edges length in the shortest path
129 : // ppi_output is pointer to a GrB_Vector, where the i-th entry is pi(i), the
130 : // parent of i-th vertex in the shortest path
131 : // ph_output is pointer to a GrB_Vector, where the i-th entry is h(s,i), the
132 : // number of edges from s to i in the shortest path
133 : // A has zeros on diagonal and weights on corresponding entries of edges
134 : // s is given index for source vertex
135 5 : GrB_Info LAGraph_BF_full
136 : (
137 : GrB_Vector *pd_output, //the pointer to the vector of distance
138 : GrB_Vector *ppi_output, //the pointer to the vector of parent
139 : GrB_Vector *ph_output, //the pointer to the vector of hops
140 : const GrB_Matrix A, //matrix for the graph
141 : const GrB_Index s //given index of the source
142 : )
143 : {
144 : GrB_Info info;
145 5 : char *msg = NULL ;
146 : // tmp vector to store distance vector after n (i.e., V) loops
147 5 : GrB_Vector d = NULL, dtmp = NULL;
148 5 : GrB_Matrix Atmp = NULL;
149 : GrB_Type BF_Tuple3;
150 :
151 : GrB_BinaryOp BF_lMIN_Tuple3;
152 : GrB_BinaryOp BF_PLUSrhs_Tuple3;
153 : GrB_BinaryOp BF_EQ_Tuple3;
154 :
155 : GrB_Monoid BF_lMIN_Tuple3_Monoid;
156 : GrB_Semiring BF_lMIN_PLUSrhs_Tuple3;
157 :
158 : GrB_Index nrows, ncols, n, nz; // n = # of row/col, nz = # of nnz in graph
159 5 : GrB_Index *I = NULL, *J = NULL; // for col/row indices of entries from A
160 5 : GrB_Index *h = NULL, *pi = NULL;
161 5 : double *w = NULL;
162 5 : BF_Tuple3_struct *W = NULL;
163 :
164 5 : LG_ASSERT (A != NULL && pd_output != NULL &&
165 : ppi_output != NULL && ph_output != NULL, GrB_NULL_POINTER) ;
166 :
167 5 : *pd_output = NULL;
168 5 : *ppi_output = NULL;
169 5 : *ph_output = NULL;
170 5 : GRB_TRY (GrB_Matrix_nrows (&nrows, A)) ;
171 5 : GRB_TRY (GrB_Matrix_ncols (&ncols, A)) ;
172 5 : GRB_TRY (GrB_Matrix_nvals (&nz, A));
173 5 : LG_ASSERT_MSG (nrows == ncols, -1002, "A must be square") ;
174 5 : n = nrows;
175 5 : LG_ASSERT_MSG (s < n, GrB_INVALID_INDEX, "invalid source node") ;
176 :
177 : //--------------------------------------------------------------------------
178 : // create all GrB_Type GrB_BinaryOp GrB_Monoid and GrB_Semiring
179 : //--------------------------------------------------------------------------
180 : // GrB_Type
181 5 : GRB_TRY (GrB_Type_new(&BF_Tuple3, sizeof(BF_Tuple3_struct)));
182 :
183 : // GrB_BinaryOp
184 5 : GRB_TRY (GrB_BinaryOp_new(&BF_EQ_Tuple3,
185 : (LAGraph_binary_function) (&BF_EQ), GrB_BOOL, BF_Tuple3, BF_Tuple3));
186 5 : GRB_TRY (GrB_BinaryOp_new(&BF_lMIN_Tuple3,
187 : (LAGraph_binary_function) (&BF_lMIN), BF_Tuple3, BF_Tuple3, BF_Tuple3));
188 5 : GRB_TRY (GrB_BinaryOp_new(&BF_PLUSrhs_Tuple3,
189 : (LAGraph_binary_function)(&BF_PLUSrhs),
190 : BF_Tuple3, BF_Tuple3, BF_Tuple3));
191 :
192 : // GrB_Monoid
193 5 : BF_Tuple3_struct BF_identity = (BF_Tuple3_struct) { .w = INFINITY,
194 : .h = UINT64_MAX, .pi = UINT64_MAX };
195 5 : GRB_TRY (GrB_Monoid_new_UDT(&BF_lMIN_Tuple3_Monoid, BF_lMIN_Tuple3,
196 : &BF_identity));
197 :
198 : //GrB_Semiring
199 5 : GRB_TRY (GrB_Semiring_new(&BF_lMIN_PLUSrhs_Tuple3,
200 : BF_lMIN_Tuple3_Monoid, BF_PLUSrhs_Tuple3));
201 :
202 : //--------------------------------------------------------------------------
203 : // allocate arrays used for tuplets
204 : //--------------------------------------------------------------------------
205 :
206 5 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &I, nz, sizeof(GrB_Index), msg)) ;
207 5 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &J, nz, sizeof(GrB_Index), msg)) ;
208 5 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &w, nz, sizeof(double), msg)) ;
209 5 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &W, nz, sizeof(BF_Tuple3_struct),
210 : msg)) ;
211 :
212 : //--------------------------------------------------------------------------
213 : // create matrix Atmp based on A, while its entries become BF_Tuple3 type
214 : //--------------------------------------------------------------------------
215 :
216 5 : GRB_TRY (GrB_Matrix_extractTuples_FP64(I, J, w, &nz, A));
217 : int nthreads, nthreads_outer, nthreads_inner ;
218 5 : LG_TRY (LAGraph_GetNumThreads (&nthreads_outer, &nthreads_inner, msg)) ;
219 5 : nthreads = nthreads_outer * nthreads_inner ;
220 5 : printf ("nthreads %d\n", nthreads) ;
221 : int64_t k;
222 : #pragma omp parallel for num_threads(nthreads) schedule(static)
223 951 : for (k = 0; k < nz; k++)
224 : {
225 946 : if (w[k] == 0) //diagonal entries
226 : {
227 184 : W[k] = (BF_Tuple3_struct) { .w = 0, .h = 0, .pi = 0 };
228 : }
229 : else
230 : {
231 762 : W[k] = (BF_Tuple3_struct) { .w = w[k], .h = 1, .pi = I[k] + 1 };
232 : }
233 : }
234 5 : GRB_TRY (GrB_Matrix_new(&Atmp, BF_Tuple3, n, n));
235 5 : GRB_TRY (GrB_Matrix_build_UDT(Atmp, I, J, W, nz, BF_lMIN_Tuple3));
236 :
237 : //--------------------------------------------------------------------------
238 : // create and initialize "distance" vector d
239 : //--------------------------------------------------------------------------
240 5 : GRB_TRY (GrB_Vector_new(&d, BF_Tuple3, n));
241 : // initial distance from s to itself
242 5 : BF_Tuple3_struct d0 = (BF_Tuple3_struct) { .w = 0, .h = 0, .pi = 0 };
243 5 : GRB_TRY (GrB_Vector_setElement_UDT(d, &d0, s));
244 :
245 : //--------------------------------------------------------------------------
246 : // start the Bellman Ford process
247 : //--------------------------------------------------------------------------
248 : // copy d to dtmp in order to create a same size of vector
249 5 : GRB_TRY (GrB_Vector_dup(&dtmp, d));
250 5 : bool same= false; // variable indicating if d == dtmp
251 5 : int64_t iter = 0; // number of iterations
252 :
253 : // terminate when no new path is found or more than V-1 loops
254 93 : while (!same && iter < n - 1)
255 : {
256 : // execute semiring on d and A, and save the result to dtmp
257 88 : GRB_TRY (GrB_vxm(dtmp, GrB_NULL, GrB_NULL, BF_lMIN_PLUSrhs_Tuple3,
258 : d, Atmp, GrB_NULL));
259 :
260 88 : LG_TRY (LAGraph_Vector_IsEqualOp (&same, dtmp, d, BF_EQ_Tuple3, NULL));
261 88 : if (!same)
262 : {
263 85 : GrB_Vector ttmp = dtmp;
264 85 : dtmp = d;
265 85 : d = ttmp;
266 : }
267 88 : iter ++;
268 : }
269 :
270 : // check for negative-weight cycle only when there was a new path in the
271 : // last loop, otherwise, there can't be a negative-weight cycle.
272 5 : if (!same)
273 : {
274 : // execute semiring again to check for negative-weight cycle
275 2 : GRB_TRY (GrB_vxm(dtmp, GrB_NULL, GrB_NULL, BF_lMIN_PLUSrhs_Tuple3,
276 : d, Atmp, GrB_NULL));
277 :
278 : // if d != dtmp, then there is a negative-weight cycle in the graph
279 2 : LG_TRY (LAGraph_Vector_IsEqualOp (&same, dtmp, d, BF_EQ_Tuple3, NULL));
280 2 : if (!same)
281 : {
282 : // printf("A negative-weight cycle found. \n");
283 2 : LG_FREE_ALL;
284 2 : return (GrB_NO_VALUE) ;
285 : }
286 : }
287 :
288 : //--------------------------------------------------------------------------
289 : // extract tuple from "distance" vector d and create GrB_Vectors for output
290 : //--------------------------------------------------------------------------
291 :
292 3 : GRB_TRY (GrB_Vector_extractTuples_UDT (I, (void *) W, &nz, d));
293 3 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &h , nz, sizeof(GrB_Index), msg)) ;
294 3 : LAGRAPH_TRY (LAGraph_Malloc ((void **) &pi, nz, sizeof(GrB_Index), msg)) ;
295 :
296 111 : for (k = 0; k < nz; k++)
297 : {
298 108 : w [k] = W[k].w ;
299 108 : h [k] = W[k].h ;
300 108 : pi[k] = W[k].pi;
301 : }
302 3 : GRB_TRY (GrB_Vector_new(pd_output, GrB_FP64, n));
303 3 : GRB_TRY (GrB_Vector_new(ppi_output, GrB_UINT64, n));
304 3 : GRB_TRY (GrB_Vector_new(ph_output, GrB_UINT64, n));
305 3 : GRB_TRY (GrB_Vector_build_FP64 (*pd_output , I, w , nz,GrB_MIN_FP64 ));
306 3 : GRB_TRY (GrB_Vector_build_UINT64(*ppi_output, I, pi, nz,GrB_MIN_UINT64));
307 3 : GRB_TRY (GrB_Vector_build_UINT64(*ph_output , I, h , nz,GrB_MIN_UINT64));
308 3 : LG_FREE_ALL;
309 3 : return (GrB_SUCCESS) ;
310 : }
|
