Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LG_check_rcc: A hand coded RCC algorithm for CSR matircies.
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 Gabriel Gomez, Texas A&M University
15 :
16 : //------------------------------------------------------------------------------
17 :
18 : #include <stdlib.h>
19 : #include "LG_internal.h"
20 : #include "LG_test.h"
21 : #include "LG_Xtest.h"
22 :
23 : #undef LG_FREE_WORK
24 : #undef LG_FREE_ALL
25 :
26 : #define LG_FREE_WORK \
27 : { \
28 : /* free any workspace used here */ \
29 : LAGraph_Free(&a_space, NULL) ; \
30 : GrB_free (&cont) ; \
31 : LAGraph_Free((void **)&Ai, NULL) ; \
32 : LAGraph_Free((void **)&Ap, NULL) ; \
33 : LAGraph_Free((void **)&slice, NULL) ; \
34 : }
35 :
36 : #define LG_FREE_ALL \
37 : { \
38 : /* free any workspace used here */ \
39 : LG_FREE_WORK ; \
40 : /* free all the output variable(s) */ \
41 : LAGraph_Free((void **)&rcc, NULL) ; \
42 : GrB_free (rccs) ; \
43 : }
44 : #define TIMINGS
45 : #ifdef TIMINGS
46 7 : static void print_timings (const double timings [16])
47 : {
48 7 : double total = timings [0] + timings [1] + timings [2] + timings [3] + timings [4];
49 7 : printf ("RCC %12.6f (%4.1f%%) init\n", timings [0], 100. * timings [0] / total) ;
50 7 : printf ("RCC %12.6f (%4.1f%%) counting edges\n", timings [1], 100. * timings [1] / total) ;
51 7 : printf ("RCC %12.6f (%4.1f%%) counting nodes\n", timings [2], 100. * timings [2] / total) ;
52 7 : printf ("RCC %12.6f (%4.1f%%) cumulative sum\n", timings [3], 100. * timings [3] / total) ;
53 7 : printf ("RCC %12.6f (%4.1f%%) calculation\n", timings [4], 100. * timings [4] / total) ;
54 7 : }
55 : #endif
56 :
57 : //Scuffed upperbound function
58 7 : static int64_t LG_binary_search // returns upperbound - 1
59 : (
60 : const int64_t pivot,
61 : const int64_t *LG_RESTRICT X_0, // search in X [p_start..p_end_-1]
62 : const int64_t p_start,
63 : const int64_t p_end
64 : )
65 : {
66 :
67 : //--------------------------------------------------------------------------
68 : // find where the Pivot appears in X
69 : //--------------------------------------------------------------------------
70 :
71 : // binary search of X [p_start...p_end-1] for the Pivot
72 7 : int64_t pleft = p_start ;
73 7 : int64_t pright = p_end;
74 76 : while (pleft < pright)
75 : {
76 69 : int64_t pmiddle = pleft + (pright - pleft) / 2 ;
77 69 : bool less = (X_0 [pmiddle] < pivot) ;
78 69 : pleft = less ? pmiddle + 1 : pleft ;
79 69 : pright = less ? pright : pmiddle ;
80 : }
81 7 : if(X_0[pleft] <= pivot)
82 7 : pleft++;
83 7 : return (--pleft) ;
84 : }
85 :
86 :
87 7 : int LG_check_rcc
88 :
89 : (
90 : // output:
91 : //rccs(i): rich club coefficent of i
92 : GrB_Vector *rccs,
93 :
94 : // input:
95 : LAGraph_Graph G, //input graph
96 : char *msg
97 : )
98 : {
99 : #if LG_SUITESPARSE_GRAPHBLAS_V10
100 7 : GxB_Container cont = NULL;
101 7 : GrB_Matrix A = G->A;
102 7 : int64_t *Ap = NULL, *Ai = NULL;
103 7 : void *a_space = NULL;
104 7 : GrB_Type p_type = NULL, i_type = NULL;
105 7 : int p_hand = 0, i_hand = 0;
106 7 : int n_threads = LG_nthreads_outer * LG_nthreads_inner;
107 7 : uint64_t p_n = 0, i_n = 0, p_size = 0, i_size = 0 ;
108 7 : int64_t max_deg = 0;
109 7 : uint64_t *epd = NULL, *vpd = NULL;
110 7 : int64_t *LG_RESTRICT slice = NULL;
111 7 : double *rcc = NULL;
112 : #ifdef TIMINGS
113 : double timings [16] ;
114 7 : memset(timings, 0, 16*sizeof(double)) ;
115 7 : double tic = LAGraph_WallClockTime ( ) ;
116 7 : LG_SET_BURBLE (false) ;
117 : #endif
118 :
119 : //--------------------------------------------------------------------------
120 : // Check inputs
121 : //--------------------------------------------------------------------------
122 7 : LG_TRY (LAGraph_CheckGraph (G, msg)) ;
123 7 : LG_ASSERT (rccs != NULL, GrB_NULL_POINTER);
124 :
125 7 : LG_ASSERT_MSG(
126 : G->kind == LAGraph_ADJACENCY_UNDIRECTED, GrB_INVALID_VALUE,
127 : "G->A must be symmetric") ;
128 7 : LG_ASSERT_MSG(
129 : G->is_symmetric_structure == LAGraph_TRUE, GrB_INVALID_VALUE,
130 : "G->A must be symmetric") ;
131 7 : LG_ASSERT_MSG (G->out_degree != NULL, GrB_EMPTY_OBJECT,
132 : "G->out_degree must be defined") ;
133 7 : LG_ASSERT_MSG (G->nself_edges == 0, GrB_INVALID_VALUE,
134 : "G->nself_edges must be zero") ;
135 7 : GRB_TRY(GxB_Container_new(&cont)) ;
136 7 : GRB_TRY(GxB_unload_Matrix_into_Container(A, cont, NULL)) ;
137 7 : LG_ASSERT_MSG(cont->format == GxB_SPARSE, GrB_NOT_IMPLEMENTED,
138 : "Matrix must be sparse") ;
139 7 : LG_TRY (LAGraph_Malloc(
140 : (void **) &Ap, cont->nrows + 1, sizeof(uint64_t), NULL)) ;
141 7 : LG_TRY (LAGraph_Malloc(
142 : (void **) &Ai, cont->nvals, sizeof(uint64_t), NULL)) ;
143 7 : p_n = cont->nrows + 1; i_n = cont->nvals;
144 7 : GRB_TRY (GrB_Vector_extractTuples_INT64(
145 : NULL, Ap, &p_n, cont->p)) ;
146 7 : GRB_TRY (GrB_Vector_extractTuples_INT64(
147 : NULL, Ai, &i_n, cont->i)) ;
148 7 : GRB_TRY (GxB_load_Matrix_from_Container(A, cont, NULL)) ;
149 7 : GRB_TRY (GrB_Vector_reduce_INT64(
150 : &max_deg, NULL, GrB_MAX_MONOID_INT64, G->out_degree, NULL)) ;
151 7 : int64_t i = 0;
152 : #ifdef TIMINGS
153 7 : timings[0] = LAGraph_WallClockTime ( );
154 : #endif
155 7 : LG_TRY (LAGraph_Calloc(&a_space, max_deg * 2, sizeof(uint64_t), NULL)) ;
156 7 : LG_TRY (LAGraph_Malloc((void **)&slice, n_threads + 1, sizeof(int64_t), NULL)) ;
157 7 : epd = a_space ;
158 7 : vpd = ((uint64_t *) a_space) + max_deg ;
159 7 : LG_TRY (LAGraph_Malloc((void **) &rcc, max_deg, sizeof(double), NULL)) ;
160 7 : LG_eslice (slice, i_n, n_threads) ;
161 : int tid;
162 : #pragma omp parallel for num_threads(n_threads) schedule(static, 1) private(i)
163 14 : for (tid = 0 ; tid < n_threads ; tid++)
164 : {
165 7 : int64_t loc_sum = 0, dp = 0;
166 : int64_t loc_arr[1024];
167 7 : memset(loc_arr, 0, 1024 * sizeof(int64_t));
168 7 : i = slice[tid];
169 7 : int64_t ptr = LG_binary_search(i, Ap, 0, p_n - 1) ;
170 8683 : while(i < slice[tid + 1])
171 : {
172 18422 : while(Ap[ptr + 1] <= i) ++ptr;
173 8676 : int64_t dp = Ap[ptr + 1] - Ap[ptr];
174 8676 : if(dp <= 1024)
175 205376 : for(; i < slice[tid + 1] && i < Ap[ptr + 1]; ++i)
176 : {
177 196701 : uint64_t di = Ap[Ai[i] + 1] - Ap[Ai[i]];
178 196701 : loc_arr[dp - 1] += (dp < di) + (dp <= di);
179 : }
180 : else
181 : {
182 1 : loc_sum = 0;
183 1830 : for(; i < slice[tid + 1] && i < Ap[ptr + 1]; ++i)
184 : {
185 1829 : uint64_t di = Ap[Ai[i]+1] - Ap[Ai[i]];
186 1829 : loc_sum += (dp < di) + (dp <= di);
187 : }
188 : #pragma omp atomic
189 1 : epd[dp - 1] += loc_sum ;
190 : }
191 : }
192 : #pragma omp critical
193 : {
194 1199 : for(int64_t j = 0; j < 1024 && j < max_deg; ++j)
195 : {
196 1192 : epd[j] += loc_arr[j];
197 : }
198 : }
199 : }
200 : #ifdef TIMINGS
201 7 : timings[1] = LAGraph_WallClockTime ( );
202 : #endif
203 :
204 : #pragma omp parallel
205 : {
206 : int64_t loc_arr[1024];
207 7 : memset(loc_arr, 0, 1024 * sizeof(int64_t));
208 : #pragma omp for schedule(static)
209 9766 : for(i = 0; i < p_n - 1; ++i)
210 : {
211 9759 : int64_t dp = Ap[i + 1] - Ap[i] - 1;
212 9759 : if(dp < 0) continue;
213 8676 : if(dp < 1024)
214 : {
215 8675 : ++loc_arr[dp];
216 : }
217 : else
218 : {
219 : #pragma omp atomic
220 1 : ++vpd[dp];
221 : }
222 : }
223 : #pragma omp critical
224 : {
225 1199 : for(int64_t j = 0; j < 1024 && j < max_deg; ++j)
226 : {
227 1192 : vpd[j] += loc_arr[j];
228 : }
229 : }
230 : }
231 :
232 : #ifdef TIMINGS
233 7 : timings[2] = LAGraph_WallClockTime ( );
234 : #endif
235 : //run a cummulative sum (backwards)
236 1997 : for(i = max_deg - 1; i > 0; --i)
237 : {
238 1990 : vpd[i-1] += vpd[i] ;
239 1990 : epd[i-1] += epd[i] ;
240 : }
241 : #ifdef TIMINGS
242 7 : timings[3] = LAGraph_WallClockTime ( );
243 : #endif
244 : #pragma omp parallel for schedule(static)
245 2004 : for(i = 0; i < max_deg; ++i)
246 : {
247 1997 : rcc[i] = ((double)epd[i]) / ((double)vpd[i] * ((double) vpd[i] - 1.0)) ;
248 : }
249 : #ifdef TIMINGS
250 7 : timings[4] = LAGraph_WallClockTime ( );
251 7 : timings[4] -= timings[3];
252 7 : timings[3] -= timings[2];
253 7 : timings[2] -= timings[1];
254 7 : timings[1] -= timings[0];
255 7 : timings[0] -= tic;
256 :
257 7 : print_timings(timings);
258 7 : LG_SET_BURBLE(false);
259 : #endif
260 7 : epd = vpd = NULL;
261 7 : GRB_TRY (GrB_Vector_new(rccs, GrB_FP64, max_deg));
262 7 : GRB_TRY (GxB_Vector_load(
263 : *rccs, (void **) &rcc, GrB_FP64, max_deg, max_deg * sizeof(double),
264 : GrB_DEFAULT, NULL)) ;
265 7 : LG_FREE_WORK ;
266 7 : return (GrB_SUCCESS) ;
267 : #else
268 : printf ("RCC not implemented for GraphBLAS versions under 10\n") ;
269 : return (GrB_NOT_IMPLEMENTED) ;
270 : #endif
271 : }
272 : #undef TIMINGS
|
