Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LAGraph_BF_basic_pushpull: Bellman-Ford method for single source shortest
3 : // paths
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 just the path lengths,
20 : // and not the parent vector.
21 :
22 : // LAGraph_BF_basic_pushpull: Bellman-Ford single source shortest paths,
23 : // returning just the shortest path lengths.
24 :
25 : // LAGraph_BF_basic_pushpull performs a Bellman-Ford to find out shortest path
26 : // length from given source vertex s in the range of [0, n) on graph given as
27 : // matrix A with size n by n. The sparse matrix A has entry A(i, j) if there is
28 : // edge from vertex i to vertex j with weight w, then A(i, j) = w. Furthermore,
29 : // LAGraph_BF_basic requires A(i, i) = 0 for all 0 <= i < n.
30 :
31 : // LAGraph_BF_basic returns GrB_SUCCESS if successful, or GrB_NO_VALUE if it
32 : // detects a negative-weight cycle. The GrB_Vector d(k) (i.e., *pd_output)
33 : // will be NULL when negative-weight cycle detected. Otherwise, the vector d
34 : // has d(k) as the shortest distance from s to k.
35 :
36 : // todo: this is a candidate for inclusion as a src/algorithm
37 :
38 : //------------------------------------------------------------------------------
39 :
40 : #define LG_FREE_ALL \
41 : { \
42 : GrB_free(&d) ; \
43 : GrB_free(&dtmp) ; \
44 : }
45 :
46 : #include <LAGraph.h>
47 : #include <LAGraphX.h>
48 : #include <LG_internal.h> // from src/utility
49 :
50 :
51 : // Given a n-by-n adjacency matrix A and a source vertex s.
52 : // If there is no negative-weight cycle reachable from s, return the distances
53 : // of shortest paths from s as vector d. Otherwise, return d=NULL if there is
54 : // negative-weight cycle.
55 : // pd_output = &d, where d is a GrB_Vector with d(k) as the shortest distance
56 : // from s to k when no negative-weight cycle detected, otherwise, d = NULL.
57 :
58 : // A must have explicit zeros on the diagonal and weights on corresponding
59 : // entries of edges. s is given index for source vertex.
60 :
61 15 : GrB_Info LAGraph_BF_basic_pushpull
62 : (
63 : GrB_Vector *pd_output, //the pointer to the vector of distance
64 : const GrB_Matrix A, //matrix for the graph
65 : const GrB_Matrix AT, //transpose of A (optional)
66 : const GrB_Index s //given index of the source
67 : )
68 : {
69 : GrB_Info info;
70 15 : char *msg = NULL ;
71 : GrB_Index nrows, ncols, nvalA;
72 : // tmp vector to store distance vector after n (i.e., V) loops
73 15 : GrB_Vector d = NULL, dtmp = NULL;
74 :
75 15 : LG_ASSERT ((A != NULL || AT != NULL) && pd_output != NULL,
76 : GrB_NULL_POINTER);
77 :
78 15 : (*pd_output) = NULL;
79 : bool use_vxm_with_A;
80 15 : if (A == NULL)
81 : {
82 5 : GRB_TRY (GrB_Matrix_nrows (&nrows, AT)) ;
83 5 : GRB_TRY (GrB_Matrix_ncols (&ncols, AT)) ;
84 5 : GRB_TRY (GrB_Matrix_nvals (&nvalA, AT)) ;
85 5 : use_vxm_with_A = false;
86 : }
87 : else
88 : {
89 10 : GRB_TRY (GrB_Matrix_nrows (&nrows, A)) ;
90 10 : GRB_TRY (GrB_Matrix_ncols (&ncols, A)) ;
91 10 : GRB_TRY (GrB_Matrix_nvals (&nvalA, A)) ;
92 10 : use_vxm_with_A = true;
93 : }
94 :
95 : // push/pull requires both A and AT
96 15 : bool push_pull = (A != NULL && AT != NULL) ;
97 :
98 15 : LG_ASSERT_MSG (nrows == ncols, -1002, "A must be square") ;
99 :
100 15 : GrB_Index n = nrows; // n = # of vertices in graph
101 : // average node degree
102 15 : double dA = (n == 0) ? 0 : (((double) nvalA) / (double) n) ;
103 :
104 15 : LG_ASSERT_MSG (s < n, GrB_INVALID_INDEX, "invalid source node") ;
105 :
106 : // values used to determine if d should be converted to dense
107 : // dthreshold is used when only A or AT is available
108 : int64_t dthreshold ;
109 15 : if (A == NULL)
110 : {
111 5 : dthreshold = LAGRAPH_MAX (256, sqrt ((double) n)) ;
112 : }
113 : else
114 : {
115 10 : dthreshold = n/2;
116 : }
117 : // refer to GraphBLAS/Source/GB_AxB_select.c
118 : // convert d to dense when GB_AxB_select intend to use Gustavson
119 15 : size_t csize = sizeof(double) + sizeof(int64_t) ;
120 15 : double gs_memory = ((double) n * (double) csize)/1e9 ;
121 :
122 15 : bool dsparse = true;
123 :
124 : // Initialize distance vector, change the d[s] to 0
125 15 : GRB_TRY (GrB_Vector_new(&d, GrB_FP64, n));
126 15 : GRB_TRY (GrB_Vector_setElement_FP64(d, 0, s));
127 : // copy d to dtmp in order to create a same size of vector
128 15 : GRB_TRY (GrB_Vector_dup(&dtmp, d));
129 :
130 15 : int64_t iter = 0; //number of iterations
131 15 : bool same = false; //variable indicating if d=dtmp
132 :
133 : // terminate when no new path is found or more than n-1 loops
134 279 : while (!same && iter < n - 1)
135 : {
136 :
137 : // double t2 = LAGraph_WallClockTime ( ) ;
138 :
139 : // excute semiring on d and A, and save the result to d
140 264 : if (!use_vxm_with_A)
141 : {
142 171 : GRB_TRY (GrB_mxv (dtmp, NULL, NULL, GrB_MIN_PLUS_SEMIRING_FP64,
143 : AT, d, NULL));
144 : }
145 : else
146 : {
147 93 : GRB_TRY (GrB_vxm (dtmp, NULL, NULL, GrB_MIN_PLUS_SEMIRING_FP64,
148 : d, A, NULL));
149 : }
150 264 : LG_TRY (LAGraph_Vector_IsEqual (&same, dtmp, d, NULL));
151 264 : if (!same)
152 : {
153 255 : GrB_Vector ttmp = dtmp;
154 255 : dtmp = d;
155 255 : d = ttmp;
156 : }
157 264 : iter++;
158 :
159 : // t2 = LAGraph_WallClockTime ( ) - t2 ;
160 :
161 : GrB_Index dnz ;
162 264 : GRB_TRY (GrB_Vector_nvals (&dnz, d)) ;
163 :
164 264 : if (dsparse)
165 : {
166 105 : if (!push_pull)
167 : {
168 100 : if (dnz > dthreshold)
169 : {
170 5 : dsparse = false;
171 : }
172 : }
173 : else
174 : {
175 5 : double heap_memory = (((double) dnz+1) *
176 5 : 5 * (double) (sizeof(int64_t))) / 1e9;
177 5 : int log2dnz = 0 ;
178 20 : while (dnz > 0)
179 : {
180 15 : dnz = dnz / 2 ;
181 15 : log2dnz++ ;
182 : }
183 5 : dsparse = (4 * log2dnz * heap_memory < gs_memory);
184 5 : use_vxm_with_A = dsparse;
185 : }
186 :
187 105 : if (!dsparse)
188 : {
189 10 : GRB_TRY (GrB_Vector_setElement_FP64(d, 1e-16, s));
190 10 : GRB_TRY (GrB_assign (d, d, NULL, INFINITY, GrB_ALL, n,
191 : GrB_DESC_C)) ;
192 10 : GRB_TRY (GrB_Vector_setElement_FP64(d, 0, s));
193 : }
194 : }
195 : }
196 :
197 : // check for negative-weight cycle only when there was a new path in the
198 : // last loop, otherwise, there can't be a negative-weight cycle.
199 15 : if (!same)
200 : {
201 : // excute semiring again to check for negative-weight cycle
202 6 : if (!use_vxm_with_A)
203 : {
204 4 : GRB_TRY (GrB_mxv(dtmp, NULL, NULL,
205 : GrB_MIN_PLUS_SEMIRING_FP64, AT, d, NULL));
206 : }
207 : else
208 : {
209 2 : GRB_TRY (GrB_vxm(dtmp, NULL, NULL,
210 : GrB_MIN_PLUS_SEMIRING_FP64, d, A, NULL));
211 : }
212 6 : LG_TRY (LAGraph_Vector_IsEqual (&same, dtmp, d, NULL));
213 :
214 : // if d != dtmp, then there is a negative-weight cycle in the graph
215 6 : if (!same)
216 : {
217 : // printf("A negative-weight cycle found. \n");
218 6 : LG_FREE_ALL;
219 6 : return (GrB_NO_VALUE) ;
220 : }
221 : }
222 :
223 : //--------------------------------------------------------------------------
224 : // todo: make d sparse
225 : //--------------------------------------------------------------------------
226 : /*if (!dsparse)
227 : {
228 : GRB_TRY (GrB_assign (d, d, NULL, d, GrB_ALL, n, GrB_DESC_R)) ;
229 : GRB_TRY (GrB_Vector_setElement_FP64(d, 0, s));
230 : GrB_Index dnz ;
231 : GRB_TRY (GrB_Vector_nvals (&dnz, d)) ;
232 : printf ("final nvals %.16g\n", (double) dnz) ;
233 : }*/
234 :
235 9 : (*pd_output) = d;
236 9 : d = NULL;
237 9 : LG_FREE_ALL;
238 9 : return (GrB_SUCCESS) ;
239 : }
|
