1 //------------------------------------------------------------------------------
2 // GB_hyper_prune: remove empty vectors from a hypersparse Ap, Ah list
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 // Removes empty vectors from a hypersparse list.  On input, *Ap and *Ah are
11 // assumed to be NULL.  The input arrays Ap_old and Ah_old are not modified,
12 // and thus can be shallow content from another matrix.  New hyperlists Ap and
13 // Ah are allocated, for nvec vectors, all nonempty.
14 
15 #include "GB.h"
16 
GB_hyper_prune(int64_t * restrict * p_Ap,size_t * p_Ap_size,int64_t * restrict * p_Ah,size_t * p_Ah_size,int64_t * p_nvec,const int64_t * Ap_old,const int64_t * Ah_old,const int64_t nvec_old,GB_Context Context)17 GrB_Info GB_hyper_prune
18 (
19     // output, not allocated on input:
20     int64_t *restrict *p_Ap, size_t *p_Ap_size,      // size nvec+1
21     int64_t *restrict *p_Ah, size_t *p_Ah_size,      // size nvec
22     int64_t *p_nvec,                // # of vectors, all nonempty
23     // input, not modified
24     const int64_t *Ap_old,          // size nvec_old+1
25     const int64_t *Ah_old,          // size nvec_old
26     const int64_t nvec_old,         // original number of vectors
27     GB_Context Context
28 )
29 {
30 
31     //--------------------------------------------------------------------------
32     // check inputs
33     //--------------------------------------------------------------------------
34 
35     ASSERT (p_Ap != NULL) ;
36     ASSERT (p_Ah != NULL) ;
37     ASSERT (p_nvec != NULL) ;
38     ASSERT (Ap_old != NULL) ;
39     ASSERT (Ah_old != NULL) ;
40     ASSERT (nvec_old >= 0) ;
41     (*p_Ap) = NULL ;    (*p_Ap_size) = 0 ;
42     (*p_Ah) = NULL ;    (*p_Ah_size) = 0 ;
43     (*p_nvec) = -1 ;
44 
45     int64_t *restrict W  = NULL ; size_t W_size  = 0 ;
46     int64_t *restrict Ap = NULL ; size_t Ap_size = 0 ;
47     int64_t *restrict Ah = NULL ; size_t Ah_size = 0 ;
48 
49     //--------------------------------------------------------------------------
50     // determine the # of threads to use
51     //--------------------------------------------------------------------------
52 
53     GB_GET_NTHREADS_MAX (nthreads_max, chunk, Context) ;
54     int nthreads = GB_nthreads (nvec_old, chunk, nthreads_max) ;
55 
56     //--------------------------------------------------------------------------
57     // allocate workspace
58     //--------------------------------------------------------------------------
59 
60     W = GB_MALLOC_WERK (nvec_old+1, int64_t, &W_size) ;
61     if (W == NULL)
62     {
63         // out of memory
64         return (GrB_OUT_OF_MEMORY) ;
65     }
66 
67     //--------------------------------------------------------------------------
68     // count the # of nonempty vectors
69     //--------------------------------------------------------------------------
70 
71     int64_t k ;
72     #pragma omp parallel for num_threads(nthreads) schedule(static)
73     for (k = 0 ; k < nvec_old ; k++)
74     {
75         // W [k] = 1 if the kth vector is nonempty; 0 if empty
76         W [k] = (Ap_old [k] < Ap_old [k+1]) ;
77     }
78 
79     int64_t nvec ;
80     GB_cumsum (W, nvec_old, &nvec, nthreads, Context) ;
81 
82     //--------------------------------------------------------------------------
83     // allocate the result
84     //--------------------------------------------------------------------------
85 
86     Ap = GB_MALLOC (nvec+1, int64_t, &Ap_size) ;
87     Ah = GB_MALLOC (nvec  , int64_t, &Ah_size) ;
88     if (Ap == NULL || Ah == NULL)
89     {
90         // out of memory
91         GB_FREE_WERK (&W, W_size) ;
92         GB_FREE (&Ap, Ap_size) ;
93         GB_FREE (&Ah, Ah_size) ;
94         return (GrB_OUT_OF_MEMORY) ;
95     }
96 
97     //--------------------------------------------------------------------------
98     // create the Ap and Ah result
99     //--------------------------------------------------------------------------
100 
101     #pragma omp parallel for num_threads(nthreads) schedule(static)
102     for (k = 0 ; k < nvec_old ; k++)
103     {
104         if (Ap_old [k] < Ap_old [k+1])
105         {
106             int64_t knew = W [k] ;
107             Ap [knew] = Ap_old [k] ;
108             Ah [knew] = Ah_old [k] ;
109         }
110     }
111 
112     Ap [nvec] = Ap_old [nvec_old] ;
113 
114     //--------------------------------------------------------------------------
115     // free workspace and return result
116     //--------------------------------------------------------------------------
117 
118     GB_FREE_WERK (&W, W_size) ;
119     (*p_Ap) = Ap ; (*p_Ap_size) = Ap_size ;
120     (*p_Ah) = Ah ; (*p_Ah_size) = Ah_size ;
121     (*p_nvec) = nvec ;
122     return (GrB_SUCCESS) ;
123 }
124 
125