1 //------------------------------------------------------------------------------
2 // GraphBLAS/Demo/Program/tri_demo.c: count triangles
3 //------------------------------------------------------------------------------
4
5 // SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2021, All Rights Reserved.
6 // SPDX-License-Identifier: Apache-2.0
7
8 //------------------------------------------------------------------------------
9
10 // Read a graph from a file and count the # of triangles using two methods.
11 // Usage:
12 //
13 // tri_demo < infile
14 // tri_demo 1 < infile
15 // tri_demo 0 nrows ncols ntuples method
16 // tri_demo 1 nx ny method
17 //
18 // Where infile has one line per edge in the graph; these have the form
19 //
20 // i j x
21 //
22 // where A(i,j)=x is performed by GrB_Matrix_build, to construct the matrix.
23 // The default file format is 0-based, but with "tri_demo 1 < infile" the
24 // matrix is assumed to be 1-based.
25
26 // The dimensions of A are assumed to be the largest row and column indices,
27 // plus one if the matrix is 1-based. This is done in read_matrix.c.
28 //
29 // For the second usage (tri_demo 0 ...), a random symmetric matrix is created
30 // of size nrows-by-ncols with ntuples edges (some will be duplicates so the
31 // actual number of edges will be slightly less). The method is 0 for
32 // setElement and 1 for build. The matrix will not have any self-edges, which
33 // cause the tricount method to fail.
34 //
35 // The 3rd usage (tri_demo 1 ...) creates a finite-element matrix on an
36 // nx-by-ny grid. Method is 0 to 3; refer to wathen.c for details.
37
38 // macro used by OK(...) to free workspace if an error occurs
39 #define FREE_ALL \
40 GxB_Scalar_free (&Thunk) ; \
41 GrB_Matrix_free (&C) ; \
42 GrB_Matrix_free (&A) ; \
43 GrB_Matrix_free (&L) ; \
44 GrB_Matrix_free (&U) ;
45
46 #include "graphblas_demos.h"
47
main(int argc,char ** argv)48 int main (int argc, char **argv)
49 {
50 GrB_Matrix C = NULL, A = NULL, L = NULL, U = NULL ;
51 GxB_Scalar Thunk = NULL ;
52 GrB_Info info ;
53 double tic [2], r1, r2 ;
54 OK (GrB_init (GrB_NONBLOCKING)) ;
55 int nthreads ;
56 OK (GxB_Global_Option_get (GxB_GLOBAL_NTHREADS, &nthreads)) ;
57 fprintf (stderr, "tri_demo: nthreads %d\n", nthreads) ;
58 printf ("--------------------------------------------------------------\n");
59
60 //--------------------------------------------------------------------------
61 // get a symmetric matrix with no self edges
62 //--------------------------------------------------------------------------
63
64 // get_matrix reads in a boolean matrix. It could easily be changed to
65 // read in int32 matrix instead, but this would affect the other GraphBLAS
66 // demos. So the time to typecast A = (int32) C is added to the read
67 // time, not the prep time for triangle counting.
68 simple_tic (tic) ;
69 OK (get_matrix (&C, argc, argv, true, true, true)) ;
70 GrB_Index n, nedges ;
71 OK (GrB_Matrix_nrows (&n, C)) ;
72
73 GrB_Type ctype = GrB_INT32 ;
74
75 // A = spones (C), and typecast to int32
76 OK (GrB_Matrix_new (&A, ctype, n, n)) ;
77 OK (GrB_Matrix_apply (A, NULL, NULL, GxB_ONE_INT32, C, NULL)) ;
78 double t_read = simple_toc (tic) ;
79 printf ("\ntotal time to read A matrix: %14.6f sec\n", t_read) ;
80 GrB_Matrix_free (&C) ;
81
82 OK (GxB_Scalar_new (&Thunk, GrB_INT64)) ;
83
84 // U = triu (A,1)
85 simple_tic (tic) ;
86 OK (GxB_Scalar_setElement_INT64 (Thunk, (int64_t) 1)) ;
87 OK (GrB_Matrix_new (&U, ctype, n, n)) ;
88 OK (GxB_Matrix_select (U, NULL, NULL, GxB_TRIU, A, Thunk, NULL)) ;
89 OK (GrB_Matrix_nvals (&nedges, U)) ;
90 printf ("\nn %.16g # edges %.16g\n", (double) n, (double) nedges) ;
91 double t_U = simple_toc (tic) ;
92 printf ("U=triu(A) time: %14.6f sec\n", t_U) ;
93
94 // L = tril (A,-1)
95 simple_tic (tic) ;
96 OK (GrB_Matrix_new (&L, ctype, n, n)) ;
97 OK (GxB_Scalar_setElement_INT64 (Thunk, (int64_t) (-1))) ;
98 OK (GxB_Matrix_select (L, NULL, NULL, GxB_TRIL, A, Thunk, NULL)) ;
99 double t_L = simple_toc (tic) ;
100 printf ("L=tril(A) time: %14.6f sec\n", t_L) ;
101 OK (GrB_Matrix_free (&A)) ;
102
103 int nthreads_max = 1 ;
104 #if defined ( _OPENMP )
105 nthreads_max = omp_get_max_threads ( ) ;
106 #endif
107
108 //--------------------------------------------------------------------------
109 // count the triangles via C<L> = L*U' (dot-produt)
110 //--------------------------------------------------------------------------
111
112 printf ("\n------------------------------------- dot product method:\n") ;
113
114 #define NTHREADS_MAX 2048
115 nthreads_max = MIN (nthreads_max, NTHREADS_MAX) ;
116
117 int64_t ntri2 [NTHREADS_MAX+1], nt = -1 ;
118 double t1 ;
119
120 for (int nthreads = 1 ; nthreads <= nthreads_max ; nthreads *= 2)
121 {
122 GxB_Global_Option_set (GxB_GLOBAL_NTHREADS, nthreads) ;
123
124 double t_dot [2] ;
125 OK (tricount (&(ntri2 [nthreads]), 5, NULL, NULL, L, U, t_dot)) ;
126
127 if (nthreads == 1)
128 {
129 printf ("# triangles %.16g\n", (double) ntri2 [nthreads]) ;
130 fprintf (stderr, "# triangles %.16g\n", (double) ntri2 [nthreads]) ;
131 nt = ntri2 [1] ;
132 }
133 if (ntri2 [nthreads] != nt)
134 {
135 printf ("error 5!\n") ;
136 fprintf (stderr, "error!\n") ;
137 exit (1) ;
138 }
139
140 printf ("L*U' time (dot): %14.6f sec", t_dot [0]) ;
141 if (nthreads == 1)
142 {
143 t1 = t_dot [0] ;
144 }
145 else
146 {
147 printf (" (nthreads: %d speedup %g)", nthreads, t1 / t_dot [0]) ;
148 }
149
150 printf ("\ntricount time: %14.6f sec (dot product method)\n",
151 t_dot [0] + t_dot [1]) ;
152 printf ("tri+prep time: %14.6f sec (incl time to compute L and U)\n",
153 t_dot [0] + t_dot [1] + t_U + t_L) ;
154
155 printf ("compute C time: %14.6f sec\n", t_dot [0]) ;
156 printf ("reduce (C) time: %14.6f sec\n", t_dot [1]) ;
157
158 r1 = 1e-6*nedges / (t_dot [0] + t_dot [1] + t_U + t_L) ;
159 r2 = 1e-6*nedges / (t_dot [0] + t_dot [1]) ;
160 printf ("rate %8.2f million edges/sec (incl time for U=triu(A))\n",r1);
161 printf ("rate %8.2f million edges/sec (just tricount itself)\n", r2);
162 fprintf (stderr, "GrB: C<L>=L*U' (dot) "
163 "rate %8.2f (w/ prep), %8.2f (tri)", r1, r2) ;
164 if (nthreads > 1) fprintf (stderr, " speedup: %6.2f", t1/ t_dot [0]) ;
165 fprintf (stderr, "\n") ;
166 }
167 if (nthreads_max > 1) fprintf (stderr, "\n") ;
168
169 //--------------------------------------------------------------------------
170 // method 6: C<U> = U*L' (dot)
171 //--------------------------------------------------------------------------
172
173 for (int nthreads = 1 ; nthreads <= nthreads_max ; nthreads *= 2)
174 {
175 GxB_Global_Option_set (GxB_GLOBAL_NTHREADS, nthreads) ;
176
177 double t_dot [2] ;
178 OK (tricount (&(ntri2 [nthreads]), 6, NULL, NULL, L, U, t_dot)) ;
179
180 // printf ("# triangles %.16g\n", (double) ntri2 [nthreads]) ;
181 // fprintf (stderr, "# triangles %.16g\n", (double) ntri2 [nthreads]) ;
182 if (ntri2 [nthreads] != nt)
183 {
184 printf ("error 6!\n") ;
185 fprintf (stderr, "error!\n") ;
186 exit (1) ;
187 }
188
189 printf ("L*U' time (dot): %14.6f sec", t_dot [0]) ;
190 if (nthreads == 1)
191 {
192 t1 = t_dot [0] ;
193 }
194 else
195 {
196 printf (" (nthreads: %d speedup %g)", nthreads, t1 / t_dot [0]) ;
197 }
198
199 printf ("\ntricount time: %14.6f sec (dot product method)\n",
200 t_dot [0] + t_dot [1]) ;
201 printf ("tri+prep time: %14.6f sec (incl time to compute L and U)\n",
202 t_dot [0] + t_dot [1] + t_U + t_L) ;
203
204 printf ("compute C time: %14.6f sec\n", t_dot [0]) ;
205 printf ("reduce (C) time: %14.6f sec\n", t_dot [1]) ;
206
207 r1 = 1e-6*nedges / (t_dot [0] + t_dot [1] + t_U + t_L) ;
208 r2 = 1e-6*nedges / (t_dot [0] + t_dot [1]) ;
209 printf ("rate %8.2f million edges/sec (incl time for U=triu(A))\n",r1);
210 printf ("rate %8.2f million edges/sec (just tricount itself)\n", r2);
211 fprintf (stderr, "GrB: C<U>=U*L' (dot) "
212 "rate %8.2f (w/ prep), %8.2f (tri)", r1, r2) ;
213 if (nthreads > 1) fprintf (stderr, " speedup: %6.2f", t1/ t_dot [0]) ;
214 fprintf (stderr, "\n") ;
215 }
216 if (nthreads_max > 1) fprintf (stderr, "\n") ;
217
218 //--------------------------------------------------------------------------
219 // count the triangles via C<L> = L*L (saxpy)
220 //--------------------------------------------------------------------------
221
222 printf ("\n----------------------------------- saxpy method:\n") ;
223
224 int64_t ntri1 [NTHREADS_MAX+1] ;
225
226 for (int nthreads = 1 ; nthreads <= nthreads_max ; nthreads *= 2)
227 {
228 GxB_Global_Option_set (GxB_GLOBAL_NTHREADS, nthreads) ;
229
230 double t_mark [2] = { 0, 0 } ;
231 OK (tricount (&ntri1 [nthreads], 3, NULL, NULL, L, NULL, t_mark)) ;
232 if (ntri1 [nthreads] != nt)
233 {
234 printf ("error 3!\n") ;
235 fprintf (stderr, "error!\n") ;
236 exit (1) ;
237 }
238
239 printf ("C<L>=L*L time (saxpy): %14.6f sec", t_mark [0]) ;
240 if (nthreads == 1)
241 {
242 t1 = t_mark [0] ;
243 }
244 else
245 {
246 printf (" (nthreads: %d speedup %g)", nthreads, t1 / t_mark [0]) ;
247 }
248
249 printf ("\ntricount time: %14.6f sec (saxpy method)\n",
250 t_mark [0] + t_mark [1]) ;
251 printf ("tri+prep time: %14.6f sec (incl time to compute L)\n",
252 t_mark [0] + t_mark [1] + t_L) ;
253
254 printf ("compute C time: %14.6f sec\n", t_mark [0]) ;
255 printf ("reduce (C) time: %14.6f sec\n", t_mark [1]) ;
256
257 r1 = 1e-6*((double)nedges) / (t_mark [0] + t_mark [1] + t_L) ;
258 r2 = 1e-6*((double)nedges) / (t_mark [0] + t_mark [1]) ;
259 printf ("rate %8.2f million edges/sec (incl time for L=tril(A))\n",r1);
260 printf ("rate %8.2f million edges/sec (just tricount itself)\n", r2);
261 fprintf (stderr, "GrB: C<L>=L*L (saxpy) "
262 "rate %8.2f (w/ prep), %8.2f (tri)", r1, r2) ;
263 if (nthreads > 1) fprintf (stderr, " speedup: %6.2f", t1/ t_mark [0]) ;
264 fprintf (stderr, "\n") ;
265 }
266
267 //--------------------------------------------------------------------------
268 // free workspace
269 //--------------------------------------------------------------------------
270
271 FREE_ALL ;
272 GrB_finalize ( ) ;
273 printf ("\n") ;
274 fprintf (stderr, "\n") ;
275 }
276
277