Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LAGraph/src/test/LG_check_sssp: stand-alone test for SSSP
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 Timothy A. Davis, Texas A&M University
15 :
16 : //------------------------------------------------------------------------------
17 :
18 : #include "LG_internal.h"
19 : #include "LG_test.h"
20 :
21 : // All computations are done in double precision
22 :
23 : typedef double LG_key_t ;
24 : typedef struct
25 : {
26 : int64_t name ;
27 : LG_key_t key ;
28 : }
29 : LG_Element ;
30 : #include "LG_heap.h"
31 :
32 : #undef LG_FREE_WORK
33 : #define LG_FREE_WORK \
34 : { \
35 : LAGraph_Free ((void **) &Heap, NULL) ; \
36 : LAGraph_Free ((void **) &Iheap, NULL) ; \
37 : LAGraph_Free ((void **) &distance, NULL) ; \
38 : LAGraph_Free ((void **) &parent, NULL) ; \
39 : LAGraph_Free ((void **) &path_length_in, NULL) ; \
40 : LAGraph_Free ((void **) &reachable, NULL) ; \
41 : LAGraph_Free ((void **) &reachable_in, NULL) ; \
42 : LAGraph_Free ((void **) &neighbor_weights, NULL) ; \
43 : LAGraph_Free ((void **) &neighbors, NULL) ; \
44 : GrB_free (&Row) ; \
45 : }
46 :
47 : #undef LG_FREE_ALL
48 : #define LG_FREE_ALL \
49 : { \
50 : LG_FREE_WORK ; \
51 : LAGraph_Free ((void **) &Ap, NULL) ; \
52 : LAGraph_Free ((void **) &Aj, NULL) ; \
53 : LAGraph_Free ((void **) &Ax, NULL) ; \
54 : }
55 :
56 : //------------------------------------------------------------------------------
57 : // test the results from SSSP
58 : //------------------------------------------------------------------------------
59 :
60 : // Because this method does on GxB_unpack on G->A, it should not be used in a
61 : // brutal memory test, unless the caller is prepared to reconstruct G->A
62 : // when the brutal test causes this method to return early.
63 :
64 483 : int LG_check_sssp
65 : (
66 : // input
67 : GrB_Vector Path_Length, // Path_Length(i) is the length of the
68 : // shortest path from src to node i.
69 : LAGraph_Graph G, // all edge weights must be > 0
70 : GrB_Index src,
71 : char *msg
72 : )
73 : {
74 :
75 : //--------------------------------------------------------------------------
76 : // check inputs
77 : //--------------------------------------------------------------------------
78 :
79 483 : GrB_Vector Row = NULL ;
80 483 : GrB_Index *Ap = NULL, *Aj = NULL, *neighbors = NULL ;
81 : GrB_Index Ap_size, Aj_size, Ax_size, n, ncols ;
82 483 : int64_t *queue = NULL, *Iheap = NULL, *parent = NULL ;
83 483 : LG_Element *Heap = NULL ;
84 :
85 483 : double *path_length_in = NULL, *distance = NULL, *neighbor_weights = NULL ;
86 483 : void *Ax = NULL ;
87 483 : bool *reachable = NULL, *reachable_in = NULL ;
88 :
89 483 : double tt = LAGraph_WallClockTime ( ) ;
90 483 : LG_TRY (LAGraph_CheckGraph (G, msg)) ;
91 483 : GRB_TRY (GrB_Matrix_nrows (&n, G->A)) ;
92 483 : GRB_TRY (GrB_Matrix_ncols (&ncols, G->A)) ;
93 : GrB_Type etype ;
94 : char atype_name [LAGRAPH_MAX_NAME_LEN] ;
95 483 : LG_TRY (LAGraph_Matrix_TypeName (atype_name, G->A, msg)) ;
96 483 : LG_TRY (LAGraph_TypeFromName (&etype, atype_name, msg)) ;
97 483 : int etypecode = 0 ;
98 483 : size_t etypesize = 0 ;
99 483 : double etypeinf = INFINITY ;
100 :
101 483 : if (etype == GrB_INT32)
102 : {
103 257 : etypecode = 0 ;
104 257 : etypesize = sizeof (int32_t) ;
105 257 : etypeinf = (double) INT32_MAX ;
106 : }
107 226 : else if (etype == GrB_INT64)
108 : {
109 48 : etypecode = 1 ;
110 48 : etypesize = sizeof (int64_t) ;
111 48 : etypeinf = (double) INT64_MAX ;
112 : }
113 178 : else if (etype == GrB_UINT32)
114 : {
115 12 : etypecode = 2 ;
116 12 : etypesize = sizeof (uint32_t) ;
117 12 : etypeinf = (double) UINT32_MAX ;
118 : }
119 166 : else if (etype == GrB_UINT64)
120 : {
121 12 : etypecode = 3 ;
122 12 : etypesize = sizeof (uint64_t) ;
123 12 : etypeinf = (double) UINT64_MAX ;
124 : }
125 154 : else if (etype == GrB_FP32)
126 : {
127 9 : etypecode = 4 ;
128 9 : etypesize = sizeof (float) ;
129 9 : etypeinf = INFINITY ;
130 : }
131 145 : else if (etype == GrB_FP64)
132 : {
133 145 : etypecode = 5 ;
134 145 : etypesize = sizeof (double) ;
135 145 : etypeinf = INFINITY ;
136 : }
137 :
138 483 : LG_ASSERT_MSG (etypesize > 0, GrB_NOT_IMPLEMENTED, "type not supported") ;
139 :
140 483 : bool print_timings = (n >= 2000) ;
141 :
142 : //--------------------------------------------------------------------------
143 : // get the contents of the Path_Length vector
144 : //--------------------------------------------------------------------------
145 :
146 483 : LG_TRY (LAGraph_Malloc ((void **) &path_length_in, n, sizeof (double),
147 : msg)) ;
148 483 : LG_TRY (LAGraph_Malloc ((void **) &reachable_in, n, sizeof (double), msg)) ;
149 62662 : for (int64_t i = 0 ; i < n ; i++)
150 : {
151 : double t ;
152 62179 : path_length_in [i] = INFINITY ;
153 62179 : int info = GrB_Vector_extractElement_FP64 (&t, Path_Length, i) ;
154 62179 : if (info == GrB_SUCCESS)
155 : {
156 62179 : path_length_in [i] = t ;
157 : }
158 62179 : reachable_in [i] = (path_length_in [i] < etypeinf) ;
159 : }
160 :
161 : //--------------------------------------------------------------------------
162 : // unpack the matrix in CSR form for SuiteSparse:GraphBLAS
163 : //--------------------------------------------------------------------------
164 :
165 483 : bool iso = false ;
166 : #if LAGRAPH_SUITESPARSE
167 : bool jumbled ;
168 483 : GRB_TRY (GxB_Matrix_unpack_CSR (G->A,
169 : &Ap, &Aj, (void **) &Ax, &Ap_size, &Aj_size, &Ax_size, &iso, &jumbled,
170 : NULL)) ;
171 : #endif
172 :
173 : //--------------------------------------------------------------------------
174 : // compute the SSSP of the graph, via Dijskstra's algorithm
175 : //--------------------------------------------------------------------------
176 :
177 483 : if (print_timings)
178 : {
179 20 : tt = LAGraph_WallClockTime ( ) - tt ;
180 20 : printf ("LG_check_sssp init time: %g sec\n", tt) ;
181 20 : tt = LAGraph_WallClockTime ( ) ;
182 : }
183 :
184 : // initializations
185 483 : LG_TRY (LAGraph_Malloc ((void **) &distance, n, sizeof (double), msg)) ;
186 483 : LG_TRY (LAGraph_Malloc ((void **) &reachable, n, sizeof (bool), msg)) ;
187 62662 : for (int64_t i = 0 ; i < n ; i++)
188 : {
189 62179 : distance [i] = INFINITY ;
190 62179 : reachable [i] = false ;
191 : }
192 483 : distance [src] = 0 ;
193 483 : reachable [src] = true ;
194 :
195 : #if !LAGRAPH_SUITESPARSE
196 : GRB_TRY (GrB_Vector_new (&Row, GrB_FP64, n)) ;
197 : LG_TRY (LAGraph_Malloc ((void **) &neighbors, n, sizeof (GrB_Index), msg)) ;
198 : LG_TRY (LAGraph_Malloc ((void **) &neighbor_weights, n, sizeof (double), msg)) ;
199 : LG_ASSERT (neighbors != NULL, GrB_OUT_OF_MEMORY) ;
200 : LG_ASSERT (neighbor_weights != NULL, GrB_OUT_OF_MEMORY) ;
201 : #endif
202 :
203 : // place all nodes in the heap (already in heap order)
204 483 : LG_TRY (LAGraph_Malloc ((void **) &Heap, (n+1), sizeof (LG_Element), msg)) ;
205 483 : LG_TRY (LAGraph_Malloc ((void **) &Iheap, n, sizeof (int64_t), msg)) ;
206 483 : LG_ASSERT (Heap != NULL && Iheap != NULL, GrB_OUT_OF_MEMORY) ;
207 483 : Heap [1].key = 0 ;
208 483 : Heap [1].name = src ;
209 483 : Iheap [src] = 1 ;
210 483 : int64_t p = 2 ;
211 62662 : for (int64_t i = 0 ; i < n ; i++)
212 : {
213 62179 : if (i != src)
214 : {
215 61696 : Heap [p].key = INFINITY ;
216 61696 : Heap [p].name = i ;
217 61696 : Iheap [i] = p ;
218 61696 : p++ ;
219 : }
220 : }
221 483 : int64_t nheap = n ;
222 483 : LG_ASSERT_MSG (LG_heap_check (Heap, Iheap, n, nheap) == 0, -2000,
223 : "invalid heap") ;
224 :
225 44388 : while (nheap > 0)
226 : {
227 : // extract the min element u from the top of the heap
228 44092 : LG_Element e = Heap [1] ;
229 44092 : int64_t u = e.name ;
230 :
231 44092 : double u_distance = e.key ;
232 : ASSERT (distance [u] == u_distance) ;
233 44092 : LG_heap_delete (1, Heap, Iheap, n, &nheap) ;
234 : ASSERT (Iheap [u] == 0) ;
235 44092 : reachable [u] = (u_distance < etypeinf) ;
236 :
237 44092 : if (n < 200)
238 : {
239 4946 : LG_ASSERT_MSG (LG_heap_check (Heap, Iheap, n, nheap) == 0, -2000,
240 : "invalid heap") ;
241 : }
242 :
243 44092 : if (u_distance == INFINITY)
244 : {
245 : // node u is not reachable, so no other nodes in the queue
246 : // are reachable either. All work is done.
247 187 : break ;
248 : }
249 :
250 : #if LAGRAPH_SUITESPARSE
251 : // directly access the indices of entries in A(u,:)
252 43905 : GrB_Index degree = Ap [u+1] - Ap [u] ;
253 43905 : GrB_Index *node_u_adjacency_list = Aj + Ap [u] ;
254 43905 : void *weights = ((char *) Ax) + ((iso ? 0 : Ap [u]) * etypesize) ;
255 : #else
256 : // extract the indices of entries in A(u,:)
257 : GrB_Index degree = n ;
258 : GRB_TRY (GrB_Col_extract (Row, NULL, NULL, G->A, GrB_ALL, n, u,
259 : GrB_DESC_T0)) ;
260 : GRB_TRY (GrB_Vector_extractTuples_FP64 (neighbors, neighbor_weights,
261 : °ree, Row)) ;
262 : GrB_Index *node_u_adjacency_list = neighbors ;
263 : double *weights = neighbor_weights ;
264 : #endif
265 :
266 : // traverse all entries in A(u,:)
267 614795 : for (int64_t k = 0 ; k < degree ; k++)
268 : {
269 : // consider edge (u,v) and its weight w
270 570890 : int64_t v = node_u_adjacency_list [k] ;
271 570890 : if (Iheap [v] == 0) continue ; // node v already in SSSP tree
272 : double w ;
273 : #if LAGRAPH_SUITESPARSE
274 269114 : switch (etypecode)
275 : {
276 134179 : default:
277 134179 : case 0: w = (( int32_t *) weights) [iso ? 0 : k] ; break ;
278 312 : case 1: w = (( int64_t *) weights) [iso ? 0 : k] ; break ;
279 51 : case 2: w = ((uint32_t *) weights) [iso ? 0 : k] ; break ;
280 51 : case 3: w = ((uint64_t *) weights) [iso ? 0 : k] ; break ;
281 90 : case 4: w = (( float *) weights) [iso ? 0 : k] ; break ;
282 134431 : case 5: w = (( double *) weights) [iso ? 0 : k] ; break ;
283 : }
284 : #else
285 : w = weights [iso ? 0 : k] ;
286 : #endif
287 :
288 269114 : LG_ASSERT_MSG (w > 0, -2002, "invalid graph (weights must be > 0)");
289 269114 : double new_distance = u_distance + w ;
290 269114 : if (distance [v] > new_distance)
291 : {
292 : // reduce the key of node v
293 60199 : distance [v] = new_distance ;
294 : // parent [v] = u ;
295 60199 : int64_t p = Iheap [v] ;
296 60199 : LG_ASSERT_MSG (Heap [p].name == v, -2000, "invalid heap") ;
297 60199 : LG_heap_decrease_key (p, new_distance, Heap, Iheap, n, nheap) ;
298 : }
299 : }
300 :
301 43905 : if (n < 200)
302 : {
303 4771 : LG_ASSERT_MSG (LG_heap_check (Heap, Iheap, n, nheap) == 0, -2000,
304 : "invalid heap") ;
305 : }
306 : }
307 :
308 483 : if (print_timings)
309 : {
310 20 : tt = LAGraph_WallClockTime ( ) - tt ;
311 20 : printf ("LG_check_sssp time: %g sec\n", tt) ;
312 20 : tt = LAGraph_WallClockTime ( ) ;
313 : }
314 :
315 : //--------------------------------------------------------------------------
316 : // repack the matrix in CSR form for SuiteSparse:GraphBLAS
317 : //--------------------------------------------------------------------------
318 :
319 : #if LAGRAPH_SUITESPARSE
320 483 : GRB_TRY (GxB_Matrix_pack_CSR (G->A,
321 : &Ap, &Aj, (void **) &Ax, Ap_size, Aj_size, Ax_size, iso, jumbled,
322 : NULL)) ;
323 : #endif
324 :
325 : //--------------------------------------------------------------------------
326 : // check the distance of each node
327 : //--------------------------------------------------------------------------
328 :
329 62662 : for (int64_t i = 0 ; i < n ; i++)
330 : {
331 62179 : bool ok = true ;
332 62179 : double err = 0 ;
333 62179 : if (isinf (distance [i]))
334 : {
335 18274 : ok = (path_length_in [i] == etypeinf || isinf (path_length_in [i]));
336 : }
337 : else
338 : {
339 43905 : err = fabs (path_length_in [i] - distance [i]) ;
340 43905 : double d = LAGRAPH_MAX (path_length_in [i], distance [i]) ;
341 43905 : if (err > 0) err = err / d ;
342 43905 : ok = (err < 1e-5) ;
343 : }
344 62179 : LG_ASSERT_MSG (ok, -2001, "invalid path length") ;
345 : }
346 :
347 : //--------------------------------------------------------------------------
348 : // check the reach
349 : //--------------------------------------------------------------------------
350 :
351 62662 : for (int64_t i = 0 ; i < n ; i++)
352 : {
353 62179 : bool ok = (reachable [i] == reachable_in [i]) ;
354 : #if 0
355 : printf ("reach [%ld]: %d %d\n", i, reachable [i], reachable_in [i]) ;
356 : if (!ok)
357 : {
358 : printf ("Hey! source %ld\n", src) ;
359 : GxB_print (G->A, 3) ;
360 : GxB_print (Path_Length, 3) ;
361 : for (int64_t i = 0 ; i < n ; i++)
362 : {
363 : printf ("check [%ld]: reach %d %d distance %g\n", i,
364 : reachable [i], reachable_in [i], distance [i]) ;
365 : }
366 :
367 : }
368 : #endif
369 62179 : LG_ASSERT_MSG (ok, -2001, "invalid reach") ;
370 : }
371 :
372 : //--------------------------------------------------------------------------
373 : // free workspace and return result
374 : //--------------------------------------------------------------------------
375 :
376 483 : LG_FREE_WORK ;
377 :
378 483 : if (print_timings)
379 : {
380 20 : tt = LAGraph_WallClockTime ( ) - tt ;
381 20 : printf ("LG_check_sssp check time: %g sec\n", tt) ;
382 : }
383 483 : return (GrB_SUCCESS) ;
384 : }
|
