File: GB_hyper_prune.c

package info (click to toggle)
suitesparse 1%3A7.10.1%2Bdfsg-1
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid, trixie
  • size: 254,920 kB
  • sloc: ansic: 1,134,743; cpp: 46,133; makefile: 4,875; fortran: 2,087; java: 1,826; sh: 996; ruby: 725; python: 495; asm: 371; sed: 166; awk: 44
file content (175 lines) | stat: -rw-r--r-- 6,388 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
//------------------------------------------------------------------------------
// GB_hyper_prune: prune empty vectors from a hypersparse matrix
//------------------------------------------------------------------------------

// SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2025, All Rights Reserved.
// SPDX-License-Identifier: Apache-2.0

//------------------------------------------------------------------------------

// On input, A->p and A->h may be shallow.  If modified, new arrays A->p and
// A->h are created, which are not shallow, and any existing hyper hash is
// freed.  If these arrays are not modified, and are shallow on input, then
// they remain shallow on output.  If new A->p and A->h arrays are constructed,
// the existing A->Y hyper_hash is freed.  A->p_is_32, A->j_is_32, and
// A->i_is_32 are unchanged in all cases.

#include "GB.h"

GrB_Info GB_hyper_prune
(
    GrB_Matrix A,               // matrix to prune
    GB_Werk Werk
)
{

    //--------------------------------------------------------------------------
    // check inputs
    //--------------------------------------------------------------------------

    ASSERT (A != NULL) ;
    ASSERT (GB_ZOMBIES_OK (A)) ;        // pattern not accessed
    ASSERT (GB_JUMBLED_OK (A)) ;
    ASSERT_MATRIX_OK (A, "A before hyper_prune", GB0) ;

    if (!GB_IS_HYPERSPARSE (A))
    { 
        // nothing to do
        return (GrB_SUCCESS) ;
    }

    //--------------------------------------------------------------------------
    // count # of empty vectors and check if pruning is needed
    //--------------------------------------------------------------------------

    // A->nvec_nonempty is needed to prune the hyperlist
    int64_t nvec_nonempty = GB_nvec_nonempty_update (A) ;
    if (nvec_nonempty == A->nvec)
    { 
        // nothing to prune
        return (GrB_SUCCESS) ;
    }
    #ifdef GB_DEBUG
    int64_t nvec_save = nvec_nonempty ;
    #endif

    //--------------------------------------------------------------------------
    // prune empty vectors
    //--------------------------------------------------------------------------

    GB_Ap_DECLARE (Ap_old, const) ; GB_Ap_PTR (Ap_old, A) ;
    GB_Ah_DECLARE (Ah_old, const) ; GB_Ah_PTR (Ah_old, A) ;

    GB_Ap_DECLARE (Ap_new, ) ; size_t Ap_new_size = 0 ;
    GB_Ah_DECLARE (Ah_new, ) ; size_t Ah_new_size = 0 ;

    GB_MDECL (W, , u) ; size_t W_size = 0 ;

    int64_t nvec_old = A->nvec ;

    size_t psize = (A->p_is_32) ? sizeof (uint32_t) : sizeof (uint64_t) ;
    size_t jsize = (A->j_is_32) ? sizeof (uint32_t) : sizeof (uint64_t) ;

    //--------------------------------------------------------------------------
    // determine the # of threads to use
    //--------------------------------------------------------------------------

    int nthreads_max = GB_Context_nthreads_max ( ) ;
    double chunk = GB_Context_chunk ( ) ;
    int nthreads = GB_nthreads (nvec_old, chunk, nthreads_max) ;

    //--------------------------------------------------------------------------
    // allocate workspace
    //--------------------------------------------------------------------------

    W = GB_MALLOC_MEMORY (nvec_old+1, jsize, &W_size) ;
    if (W == NULL)
    { 
        // out of memory
        return (GrB_OUT_OF_MEMORY) ;
    }
    GB_IPTR (W, A->j_is_32) ;

    //--------------------------------------------------------------------------
    // count the # of nonempty vectors and mark their locations in W
    //--------------------------------------------------------------------------

    int64_t k ;
    #pragma omp parallel for num_threads(nthreads) schedule(static)
    for (k = 0 ; k < nvec_old ; k++)
    { 
        // W [k] = 1 if the kth vector is nonempty; 0 if empty
        int nonempty = (GB_IGET (Ap_old, k) < GB_IGET (Ap_old, k+1)) ;
        // W [k] = nonempty ;
        GB_ISET (W, k, nonempty) ;
    }

    int64_t nvec_new ;
    GB_cumsum (W, A->j_is_32, nvec_old, &nvec_new, nthreads, Werk) ;

    //--------------------------------------------------------------------------
    // allocate the result
    //--------------------------------------------------------------------------

    int64_t plen_new = GB_IMAX (1, nvec_new) ;
    Ap_new = GB_MALLOC_MEMORY (plen_new+1, psize, &Ap_new_size) ;
    Ah_new = GB_MALLOC_MEMORY (plen_new  , jsize, &Ah_new_size) ;
    if (Ap_new == NULL || Ah_new == NULL)
    { 
        // out of memory
        GB_FREE_MEMORY (&W, W_size) ;
        GB_FREE_MEMORY (&Ap_new, Ap_new_size) ;
        GB_FREE_MEMORY (&Ah_new, Ah_new_size) ;
        return (GrB_OUT_OF_MEMORY) ;
    }
    GB_IPTR (Ap_new, A->p_is_32) ;
    GB_IPTR (Ah_new, A->j_is_32) ;

    //--------------------------------------------------------------------------
    // create the Ap_new and Ah_new result
    //--------------------------------------------------------------------------

    #pragma omp parallel for num_threads(nthreads) schedule(static)
    for (k = 0 ; k < nvec_old ; k++)
    {
        uint64_t p = GB_IGET (Ap_old, k) ;
        if (p < GB_IGET (Ap_old, k+1))
        { 
            uint64_t j = GB_IGET (Ah_old, k) ;
            uint64_t knew = GB_IGET (W, k) ;
            // Ap_new [knew] = p ;
            GB_ISET (Ap_new, knew, p) ;
            // Ah_new [knew] = j ;
            GB_ISET (Ah_new, knew, j) ;
        }
    }

    // Ap_new [nvec_new] = Ap_old [nvec_old] ;
    uint64_t nvals = A->nvals ;
    ASSERT (nvals == GB_IGET (Ap_old, nvec_old)) ;
    GB_ISET (Ap_new, nvec_new, nvals) ;

    //--------------------------------------------------------------------------
    // free workspace and old matrix components, including the A->Y hyper_hash
    //--------------------------------------------------------------------------

    GB_FREE_MEMORY (&W, W_size) ;
    GB_phy_free (A) ;

    //--------------------------------------------------------------------------
    // transplant the new hyperlist into A
    //--------------------------------------------------------------------------

    A->p = Ap_new ; A->p_size = Ap_new_size ;
    A->h = Ah_new ; A->h_size = Ah_new_size ;
    A->nvec = nvec_new ;
    A->plen = plen_new ;
    ASSERT (nvec_new == nvec_save) ;
    GB_nvec_nonempty_set (A, nvec_new) ;
    A->nvals = nvals ;
    A->magic = GB_MAGIC ;

    ASSERT_MATRIX_OK (A, "A after hyper_prune", GB0) ;
    return (GrB_SUCCESS) ;
}