LCOV - code coverage report
Current view: top level - src/algorithm - LAGr_TriangleCount.c (source / functions) Hit Total Coverage
Test: LAGraph code coverage report. Commit id: 7b9d2ee. Current time (UTC): 2025-06-03T21:57:17Z Lines: 103 103 100.0 %
Date: 2025-06-03 22:02:40 Functions: 2 2 100.0 %

          Line data    Source code
       1             : //------------------------------------------------------------------------------
       2             : // LAGr_TriangleCount: Triangle counting using various methods
       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             : // Count the number of triangles in a graph,
      19             : 
      20             : // This is an Advanced algorithm (G->nself_edges, G->out_degree,
      21             : // G->is_symmetric_structure are required).
      22             : 
      23             : // Given a symmetric graph A with no-self edges, LAGr_TriangleCount counts the
      24             : // number of triangles in the graph.  A triangle is a clique of size three,
      25             : // that is, 3 nodes that are all pairwise connected.
      26             : 
      27             : // One of 6 methods are used, defined below where L and U are the strictly
      28             : // lower and strictly upper triangular parts of the symmetrix matrix A,
      29             : // respectively.  Each method computes the same result, ntri:
      30             : 
      31             : //  0:  default:    use the default method (currently method Sandia_LUT)
      32             : //  1:  Burkhardt:  ntri = sum (sum ((A^2) .* A)) / 6
      33             : //  2:  Cohen:      ntri = sum (sum ((L * U) .* A)) / 2
      34             : //  3:  Sandia_LL:  ntri = sum (sum ((L * L) .* L))
      35             : //  4:  Sandia_UU:  ntri = sum (sum ((U * U) .* U))
      36             : //  5:  Sandia_LUT: ntri = sum (sum ((L * U') .* L)).  Note that L=U'.
      37             : //  6:  Sandia_ULT: ntri = sum (sum ((U * L') .* U)).  Note that U=L'.
      38             : 
      39             : // A is a square symmetric matrix, of any type.  Its values are ignored.
      40             : // Results are undefined for methods 1 and 2 if self-edges exist in A.  Results
      41             : // are undefined for all methods if A is unsymmetric.
      42             : 
      43             : // The Sandia_* methods all tend to be faster than the Burkhardt or Cohen
      44             : // methods.  For the largest graphs, Sandia_LUT tends to be fastest, except for
      45             : // the GAP-urand matrix, where the saxpy-based Sandia_LL method (L*L.*L) is
      46             : // fastest.  For many small graphs, the saxpy-based Sandia_LL and Sandia_UU
      47             : // methods are often faster that the dot-product-based methods.
      48             : 
      49             : // Reference for the Burkhardt method:  Burkhardt, Paul. "Graphing Trillions of
      50             : // Triangles." Information Visualization 16, no. 3 (July 2017): 157–66.
      51             : // https://doi.org/10.1177/1473871616666393.
      52             : 
      53             : // Reference for the Cohen method:  J. Cohen, "Graph twiddling in a mapreduce
      54             : // world," Computing in Science & Engineering, vol. 11, no. 4, pp. 29–41, 2009.
      55             : // https://doi.org/10.1109/MCSE.2009.120
      56             : 
      57             : // Reference for the "Sandia_*" methods: Wolf, Deveci, Berry, Hammond,
      58             : // Rajamanickam, "Fast linear algebra- based triangle counting with
      59             : // KokkosKernels", IEEE HPEC'17, https://dx.doi.org/10.1109/HPEC.2017.8091043
      60             : 
      61             : #define LG_FREE_ALL             \
      62             : {                               \
      63             :     GrB_free (L) ;              \
      64             :     GrB_free (U) ;              \
      65             : }
      66             : 
      67             : #include "LG_internal.h"
      68             : 
      69             : //------------------------------------------------------------------------------
      70             : // tricount_prep: construct L and U for LAGr_TriangleCount
      71             : //------------------------------------------------------------------------------
      72             : 
      73        4235 : static int tricount_prep
      74             : (
      75             :     GrB_Matrix *L,      // if present, compute L = tril (A,-1)
      76             :     GrB_Matrix *U,      // if present, compute U = triu (A, 1)
      77             :     GrB_Matrix A,       // input matrix
      78             :     char *msg
      79             : )
      80             : {
      81             :     GrB_Index n ;
      82        4235 :     GRB_TRY (GrB_Matrix_nrows (&n, A)) ;
      83             : 
      84        4235 :     if (L != NULL)
      85             :     {
      86             :         // L = tril (A,-1)
      87        3746 :         GRB_TRY (GrB_Matrix_new (L, GrB_BOOL, n, n)) ;
      88        3263 :         GRB_TRY (GrB_select (*L, NULL, NULL, GrB_TRIL, A, (int64_t) (-1),
      89             :             NULL)) ;
      90        2378 :         GRB_TRY (GrB_Matrix_wait (*L, GrB_MATERIALIZE)) ;
      91             :     }
      92             : 
      93        2867 :     if (U != NULL)
      94             :     {
      95             :         // U = triu (A,1)
      96        2626 :         GRB_TRY (GrB_Matrix_new (U, GrB_BOOL, n, n)) ;
      97        2116 :         GRB_TRY (GrB_select (*U, NULL, NULL, GrB_TRIU, A, (int64_t) 1, NULL)) ;
      98        1217 :         GRB_TRY (GrB_Matrix_wait (*U, GrB_MATERIALIZE)) ;
      99             :     }
     100        1458 :     return (GrB_SUCCESS) ;
     101             : }
     102             : 
     103             : //------------------------------------------------------------------------------
     104             : // LAGraph_tricount: count the number of triangles in a graph
     105             : //------------------------------------------------------------------------------
     106             : 
     107             : #undef  LG_FREE_ALL
     108             : #define LG_FREE_ALL                         \
     109             : {                                           \
     110             :     GrB_free (&C) ;                         \
     111             :     GrB_free (&L) ;                         \
     112             :     GrB_free (&T) ;                         \
     113             :     GrB_free (&U) ;                         \
     114             :     LAGraph_Free ((void **) &P, NULL) ;     \
     115             : }
     116             : 
     117        5856 : int LAGr_TriangleCount
     118             : (
     119             :     // output:
     120             :     uint64_t *ntriangles,
     121             :     // input:
     122             :     const LAGraph_Graph G,
     123             :     LAGr_TriangleCount_Method *p_method,
     124             :     LAGr_TriangleCount_Presort *p_presort,
     125             :     char *msg
     126             : )
     127             : {
     128             : 
     129             :     //--------------------------------------------------------------------------
     130             :     // check inputs
     131             :     //--------------------------------------------------------------------------
     132             : 
     133        5856 :     LG_CLEAR_MSG ;
     134        5856 :     GrB_Matrix C = NULL, L = NULL, U = NULL, T = NULL ;
     135        5856 :     int64_t *P = NULL ;
     136             : 
     137             :     // get the method
     138             :     LAGr_TriangleCount_Method method ;
     139        5856 :     method = (p_method == NULL) ? LAGr_TriangleCount_AutoMethod : (*p_method) ;
     140        5856 :     LG_ASSERT_MSG (
     141             :     method == LAGr_TriangleCount_AutoMethod ||  // 0: use auto method
     142             :     method == LAGr_TriangleCount_Burkhardt  ||  // 1: sum (sum ((A^2) .* A))/6
     143             :     method == LAGr_TriangleCount_Cohen      ||  // 2: sum (sum ((L * U) .*A))/2
     144             :     method == LAGr_TriangleCount_Sandia_LL  ||  // 3: sum (sum ((L * L) .* L))
     145             :     method == LAGr_TriangleCount_Sandia_UU  ||  // 4: sum (sum ((U * U) .* U))
     146             :     method == LAGr_TriangleCount_Sandia_LUT ||  // 5: sum (sum ((L * U') .* L))
     147             :     method == LAGr_TriangleCount_Sandia_ULT,    // 6: sum (sum ((U * L') .* U))
     148             :     GrB_INVALID_VALUE, "method is invalid") ;
     149             : 
     150             :     // get the presort
     151             :     LAGr_TriangleCount_Presort presort ;
     152        5846 :     presort = (p_presort == NULL) ? LAGr_TriangleCount_AutoSort : (*p_presort) ;
     153        5846 :     LG_ASSERT_MSG (
     154             :     presort == LAGr_TriangleCount_NoSort     ||
     155             :     presort == LAGr_TriangleCount_Ascending  ||
     156             :     presort == LAGr_TriangleCount_Descending ||
     157             :     presort == LAGr_TriangleCount_AutoSort,
     158             :     GrB_INVALID_VALUE, "presort is invalid") ;
     159             : 
     160        5836 :     LG_TRY (LAGraph_CheckGraph (G, msg)) ;
     161        5836 :     LG_ASSERT (ntriangles != NULL, GrB_NULL_POINTER) ;
     162        5836 :     LG_ASSERT (G->nself_edges == 0, LAGRAPH_NO_SELF_EDGES_ALLOWED) ;
     163             : 
     164        5836 :     LG_ASSERT_MSG ((G->kind == LAGraph_ADJACENCY_UNDIRECTED ||
     165             :        (G->kind == LAGraph_ADJACENCY_DIRECTED &&
     166             :         G->is_symmetric_structure == LAGraph_TRUE)),
     167             :         LAGRAPH_SYMMETRIC_STRUCTURE_REQUIRED,
     168             :         "G->A must be known to be symmetric") ;
     169             : 
     170        5836 :     if (method == LAGr_TriangleCount_AutoMethod)
     171             :     {
     172             :         // AutoMethod: use default, Sandia_LUT: sum (sum ((L * U') .* L))
     173        1433 :         method = LAGr_TriangleCount_Sandia_LUT ;
     174             :     }
     175             : 
     176             :     // only the Sandia_* methods can benefit from the presort
     177        5836 :     bool method_can_use_presort =
     178        5118 :     method == LAGr_TriangleCount_Sandia_LL || // sum (sum ((L * L) .* L))
     179        4402 :     method == LAGr_TriangleCount_Sandia_UU || // sum (sum ((U * U) .* U))
     180       10954 :     method == LAGr_TriangleCount_Sandia_LUT || // sum (sum ((L * U') .* L))
     181             :     method == LAGr_TriangleCount_Sandia_ULT ; // sum (sum ((U * L') .* U))
     182             : 
     183        5836 :     GrB_Matrix A = G->A ;
     184        5836 :     GrB_Vector Degree = G->out_degree ;
     185             : 
     186        5836 :     bool auto_sort = (presort == LAGr_TriangleCount_AutoSort) ;
     187        5836 :     if (auto_sort && method_can_use_presort)
     188             :     {
     189        1796 :         LG_ASSERT_MSG (Degree != NULL,
     190             :             LAGRAPH_NOT_CACHED, "G->out_degree is required") ;
     191             :     }
     192             : 
     193             :     //--------------------------------------------------------------------------
     194             :     // initializations
     195             :     //--------------------------------------------------------------------------
     196             : 
     197             :     GrB_Index n ;
     198        5832 :     GRB_TRY (GrB_Matrix_nrows (&n, A)) ;
     199        5832 :     GRB_TRY (GrB_Matrix_new (&C, GrB_INT64, n, n)) ;
     200        5162 :     GrB_Semiring semiring = LAGraph_plus_one_int64 ;
     201        5162 :     GrB_Monoid monoid = GrB_PLUS_MONOID_INT64 ;
     202             : 
     203             :     //--------------------------------------------------------------------------
     204             :     // heuristic sort rule
     205             :     //--------------------------------------------------------------------------
     206             : 
     207        5162 :     if (!method_can_use_presort)
     208             :     {
     209             :         // no sorting for the Burkhardt and Cohen methods: presort parameter
     210             :         // is ignored.
     211         929 :         presort = LAGr_TriangleCount_NoSort ;
     212             :     }
     213        4233 :     else if (auto_sort)
     214             :     {
     215             :         // auto selection of sorting method for Sandia_* methods
     216        1602 :         presort = LAGr_TriangleCount_NoSort ; // default is not to sort
     217             : 
     218        1602 :         if (method_can_use_presort)
     219             :         {
     220             :             // This rule is very similar to Scott Beamer's rule in the GAP TC
     221             :             // benchmark, except that it is extended to handle the ascending
     222             :             // sort needed by methods 3 and 5.  It also uses a stricter rule,
     223             :             // since the performance of triangle counting in SuiteSparse:
     224             :             // GraphBLAS is less sensitive to the sorting as compared to the
     225             :             // GAP algorithm.  This is because the dot products in SuiteSparse:
     226             :             // GraphBLAS use binary search if one vector is very sparse
     227             :             // compared to the other.  As a result, SuiteSparse:GraphBLAS needs
     228             :             // the sort for fewer matrices, as compared to the GAP algorithm.
     229             : 
     230             :             // With this rule, the GAP-kron and GAP-twitter matrices are
     231             :             // sorted, and the others remain unsorted.  With the rule in the
     232             :             // GAP tc.cc benchmark, GAP-kron and GAP-twitter are sorted, and so
     233             :             // is GAP-web, but GAP-web is not sorted here.
     234             : 
     235             :             #define NSAMPLES 2000
     236             :             GrB_Index nvals ;
     237        1614 :             GRB_TRY (GrB_Matrix_nvals (&nvals, A)) ;
     238        1602 :             if (n > NSAMPLES && ((double) nvals / ((double) n)) >= 10)
     239             :             {
     240             :                 // estimate the mean and median degrees
     241             :                 double mean, median ;
     242         479 :                 LG_TRY (LAGr_SampleDegree (&mean, &median,
     243             :                     G, true, NSAMPLES, n, msg)) ;
     244             :                 // sort if the average degree is very high vs the median
     245         467 :                 if (mean > 3 * median)
     246             :                 {
     247         292 :                     switch (method)
     248             :                     {
     249          40 :                         case LAGr_TriangleCount_Sandia_LL:
     250             :                             // 3:sum (sum ((L * L) .* L))
     251          40 :                             presort = LAGr_TriangleCount_Ascending  ;
     252          40 :                             break ;
     253          38 :                         case LAGr_TriangleCount_Sandia_UU:
     254             :                             // 4: sum (sum ((U * U) .* U))
     255          38 :                             presort = LAGr_TriangleCount_Descending ;
     256          38 :                             break ;
     257         172 :                         default:
     258             :                         case LAGr_TriangleCount_Sandia_LUT:
     259             :                             // 5: sum (sum ((L * U') .* L))
     260         172 :                             presort = LAGr_TriangleCount_Ascending  ;
     261         172 :                             break ;
     262          42 :                         case LAGr_TriangleCount_Sandia_ULT:
     263             :                             // 6: sum (sum ((U * L') .* U))
     264          42 :                             presort = LAGr_TriangleCount_Descending ;
     265          42 :                             break ;
     266             :                     }
     267             :                 }
     268             :             }
     269             :         }
     270             :     }
     271             : 
     272             :     //--------------------------------------------------------------------------
     273             :     // sort the input matrix, if requested
     274             :     //--------------------------------------------------------------------------
     275             : 
     276        5150 :     if (presort != LAGr_TriangleCount_NoSort)
     277             :     {
     278             :         // P = permutation that sorts the rows by their degree
     279        1888 :         LG_TRY (LAGr_SortByDegree (&P, G, true,
     280             :             presort == LAGr_TriangleCount_Ascending, msg)) ;
     281             : 
     282             :         // T = A (P,P) and typecast to boolean
     283        1709 :         GRB_TRY (GrB_Matrix_new (&T, GrB_BOOL, n, n)) ;
     284        1538 :         GRB_TRY (GrB_extract (T, NULL, NULL, A, (GrB_Index *) P, n,
     285             :             (GrB_Index *) P, n, NULL)) ;
     286        1184 :         A = T ;
     287             : 
     288             :         // free workspace
     289        1184 :         LG_TRY (LAGraph_Free ((void **) &P, NULL)) ;
     290             :     }
     291             : 
     292             :     //--------------------------------------------------------------------------
     293             :     // count triangles
     294             :     //--------------------------------------------------------------------------
     295             : 
     296             :     int64_t ntri ;
     297             : 
     298        4446 :     switch (method)
     299             :     {
     300             : 
     301         211 :         case LAGr_TriangleCount_Burkhardt:  // 1: sum (sum ((A^2) .* A)) / 6
     302             : 
     303         211 :             GRB_TRY (GrB_mxm (C, A, NULL, semiring, A, A, GrB_DESC_S)) ;
     304          61 :             GRB_TRY (GrB_reduce (&ntri, NULL, monoid, C, NULL)) ;
     305          61 :             ntri /= 6 ;
     306          61 :             break ;
     307             : 
     308         718 :         case LAGr_TriangleCount_Cohen: // 2: sum (sum ((L * U) .* A)) / 2
     309             : 
     310         718 :             LG_TRY (tricount_prep (&L, &U, A, msg)) ;
     311         214 :             GRB_TRY (GrB_mxm (C, A, NULL, semiring, L, U, GrB_DESC_S)) ;
     312          61 :             GRB_TRY (GrB_reduce (&ntri, NULL, monoid, C, NULL)) ;
     313          61 :             ntri /= 2 ;
     314          61 :             break ;
     315             : 
     316         489 :         case LAGr_TriangleCount_Sandia_LL: // 3: sum (sum ((L * L) .* L))
     317             : 
     318             :             // using the masked saxpy3 method
     319         489 :             LG_TRY (tricount_prep (&L, NULL, A, msg)) ;
     320         241 :             GRB_TRY (GrB_mxm (C, L, NULL, semiring, L, L, GrB_DESC_S)) ;
     321          62 :             GRB_TRY (GrB_reduce (&ntri, NULL, monoid, C, NULL)) ;
     322          62 :             break ;
     323             : 
     324         489 :         case LAGr_TriangleCount_Sandia_UU: // 4: sum (sum ((U * U) .* U))
     325             : 
     326             :             // using the masked saxpy3 method
     327         489 :             LG_TRY (tricount_prep (NULL, &U, A, msg)) ;
     328         241 :             GRB_TRY (GrB_mxm (C, U, NULL, semiring, U, U, GrB_DESC_S)) ;
     329          62 :             GRB_TRY (GrB_reduce (&ntri, NULL, monoid, C, NULL)) ;
     330          62 :             break ;
     331             : 
     332        1837 :         default:
     333             :         case LAGr_TriangleCount_Sandia_LUT: // 5: sum (sum ((L * U') .* L))
     334             : 
     335             :             // This tends to be the fastest method for most large matrices, but
     336             :             // the Sandia_ULT method is also very fast.
     337             : 
     338             :             // using the masked dot product
     339        1837 :             LG_TRY (tricount_prep (&L, &U, A, msg)) ;
     340         554 :             GRB_TRY (GrB_mxm (C, L, NULL, semiring, L, U, GrB_DESC_ST1)) ;
     341         164 :             GRB_TRY (GrB_reduce (&ntri, NULL, monoid, C, NULL)) ;
     342         164 :             break ;
     343             : 
     344         702 :         case LAGr_TriangleCount_Sandia_ULT: // 6: sum (sum ((U * L') .* U))
     345             : 
     346             :             // using the masked dot product
     347         702 :             LG_TRY (tricount_prep (&L, &U, A, msg)) ;
     348         208 :             GRB_TRY (GrB_mxm (C, U, NULL, semiring, U, L, GrB_DESC_ST1)) ;
     349          62 :             GRB_TRY (GrB_reduce (&ntri, NULL, monoid, C, NULL)) ;
     350          62 :             break ;
     351             :     }
     352             : 
     353             :     //--------------------------------------------------------------------------
     354             :     // return result
     355             :     //--------------------------------------------------------------------------
     356             : 
     357         472 :     LG_FREE_ALL ;
     358         472 :     if (p_method != NULL) (*p_method) = method ;
     359         472 :     if (p_presort != NULL) (*p_presort) = presort ;
     360         472 :     (*ntriangles) = (uint64_t) ntri ;
     361         472 :     return (GrB_SUCCESS) ;
     362             : }

Generated by: LCOV version 1.14