File: GB_convert_sparse_to_hyper.c

package info (click to toggle)
suitesparse-graphblas 7.4.0%2Bdfsg-2
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid, trixie
  • size: 67,112 kB
  • sloc: ansic: 1,072,243; cpp: 8,081; sh: 512; makefile: 506; asm: 369; python: 125; awk: 10
file content (183 lines) | stat: -rw-r--r-- 7,329 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
176
177
178
179
180
181
182
183
//------------------------------------------------------------------------------
// GB_convert_sparse_to_hyper: convert a matrix from sparse to hyperspasre
//------------------------------------------------------------------------------

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

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

// On input, the matrix may have shallow A->p content; it is safely removed.
// On output, the matrix is always hypersparse (even if out of memory).  If the
// input matrix is non-hypersparse, it is given new A->p and A->h that are not
// shallow.  If the input matrix is already hypersparse, nothing is changed
// (and in that case A->p and A->h remain shallow on output if shallow on
// input). The A->x and A->i content is not changed; it remains in whatever
// shallow/non-shallow/iso property that it had on input).

// If an out-of-memory condition occurs, all content of the matrix is cleared.

// If the input matrix A is hypersparse, bitmap or full, it is unchanged.

#include "GB.h"

GrB_Info GB_convert_sparse_to_hyper // convert from sparse to hypersparse
(
    GrB_Matrix A,           // matrix to convert to hypersparse
    GB_Context Context
)
{

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

    ASSERT_MATRIX_OK (A, "A converting to hypersparse", GB0) ;
    int64_t anz = GB_nnz (A) ;
    ASSERT (GB_ZOMBIES_OK (A)) ;
    ASSERT (GB_JUMBLED_OK (A)) ;
    ASSERT (GB_PENDING_OK (A)) ;

    //--------------------------------------------------------------------------
    // convert A from sparse to hypersparse
    //--------------------------------------------------------------------------

    if (GB_IS_SPARSE (A))
    {

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

        int64_t n = A->vdim ;
        GB_BURBLE_N (n, "(sparse to hyper) ") ;
        GB_GET_NTHREADS_MAX (nthreads_max, chunk, Context) ;
        int nthreads = GB_nthreads (n, chunk, nthreads_max) ;
        int ntasks = (nthreads == 1) ? 1 : (8 * nthreads) ;
        ntasks = GB_IMIN (ntasks, n) ;
        ntasks = GB_IMAX (ntasks, 1) ;

        //----------------------------------------------------------------------
        // count the number of non-empty vectors in A in each slice
        //----------------------------------------------------------------------

        ASSERT (A->nvec == A->plen && A->plen == n) ;

        const int64_t *restrict Ap_old = A->p ;
        size_t Ap_old_size = A->p_size ;
        bool Ap_old_shallow = A->p_shallow ;

        GB_WERK_DECLARE (Count, int64_t) ;
        GB_WERK_PUSH (Count, ntasks+1, int64_t) ;
        if (Count == NULL)
        { 
            // out of memory
            return (GrB_OUT_OF_MEMORY) ;
        }

        int tid ;
        #pragma omp parallel for num_threads(nthreads) schedule(dynamic,1)
        for (tid = 0 ; tid < ntasks ; tid++)
        {
            int64_t jstart, jend, my_nvec_nonempty = 0 ; ;
            GB_PARTITION (jstart, jend, n, tid, ntasks) ;
            for (int64_t j = jstart ; j < jend ; j++)
            { 
                if (Ap_old [j] < Ap_old [j+1]) my_nvec_nonempty++ ;
            }
            Count [tid] = my_nvec_nonempty ;
        }

        //----------------------------------------------------------------------
        // compute cumulative sum of Counts and nvec_nonempty
        //----------------------------------------------------------------------

        GB_cumsum (Count, ntasks, NULL, 1, NULL) ;
        int64_t nvec_nonempty = Count [ntasks] ;
        A->nvec_nonempty = nvec_nonempty ;

        //----------------------------------------------------------------------
        // allocate the new A->p and A->h
        //----------------------------------------------------------------------

        int64_t *restrict Ap_new = NULL ; size_t Ap_new_size = 0 ;
        int64_t *restrict Ah_new = NULL ; size_t Ah_new_size = 0 ;
        int64_t plen_new = (n == 1) ? 1 : nvec_nonempty ;
        Ap_new = GB_MALLOC (plen_new+1, int64_t, &Ap_new_size) ;
        Ah_new = GB_MALLOC (plen_new  , int64_t, &Ah_new_size) ;
        if (Ap_new == NULL || Ah_new == NULL)
        { 
            // out of memory
            GB_WERK_POP (Count, int64_t) ;
            GB_FREE (&Ap_new, Ap_new_size) ;
            GB_FREE (&Ah_new, Ah_new_size) ;
            return (GrB_OUT_OF_MEMORY) ;
        }

        //----------------------------------------------------------------------
        // transplant the new A->p and A->h into the matrix
        //----------------------------------------------------------------------

        A->plen = plen_new ;
        A->nvec = nvec_nonempty ;
        A->p = Ap_new ; A->p_size = Ap_new_size ;
        A->h = Ah_new ; A->h_size = Ah_new_size ;
        A->p_shallow = false ;
        A->h_shallow = false ;

        //----------------------------------------------------------------------
        // construct the new hyperlist in the new A->p and A->h
        //----------------------------------------------------------------------

        #pragma omp parallel for num_threads(nthreads) schedule(dynamic,1)
        for (tid = 0 ; tid < ntasks ; tid++)
        {
            int64_t jstart, jend, k = Count [tid] ;
            GB_PARTITION (jstart, jend, n, tid, ntasks) ;
            for (int64_t j = jstart ; j < jend ; j++)
            {
                if (Ap_old [j] < Ap_old [j+1])
                { 
                    // vector index j is the kth vector in the new Ah
                    Ap_new [k] = Ap_old [j] ;
                    Ah_new [k] = j ;
                    k++ ;
                }
            }
            ASSERT (k == Count [tid+1]) ;
        }

        Ap_new [nvec_nonempty] = anz ;
        A->magic = GB_MAGIC ;

        //----------------------------------------------------------------------
        // free workspace, and free the old A->p unless it's shallow
        //----------------------------------------------------------------------

        GB_WERK_POP (Count, int64_t) ;
        if (!Ap_old_shallow)
        { 
            GB_FREE (&Ap_old, Ap_old_size) ;
        }

        //----------------------------------------------------------------------
        // A is now hypersparse, but A->Y is not yet constructed
        //----------------------------------------------------------------------

        ASSERT (GB_IS_HYPERSPARSE (A)) ;
        ASSERT (A->Y == NULL && A->Y_shallow == false) ;
    }

    //--------------------------------------------------------------------------
    // A is now in hypersparse form (or left as full or bitmap)
    //--------------------------------------------------------------------------

    ASSERT (anz == GB_nnz (A)) ;
    ASSERT_MATRIX_OK (A, "A conv to hypersparse (or left full/bitmap)", GB0) ;
    ASSERT (!GB_IS_SPARSE (A)) ;
    ASSERT (GB_ZOMBIES_OK (A)) ;
    ASSERT (GB_JUMBLED_OK (A)) ;
    ASSERT (GB_PENDING_OK (A)) ;
    return (GrB_SUCCESS) ;
}