1 //------------------------------------------------------------------------------
2 // GB_convert_hyper_to_sparse: convert a matrix from hypersparse to sparse
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 // On input, the matrix may have shallow A->p and A->h content; it is safely
11 // removed. On output, the matrix is always non-hypersparse (even if out of
12 // memory). If the input matrix is hypersparse, it is given a new A->p that is
13 // not shallow. If the input matrix is already non-hypersparse, nothing is
14 // changed (and in that case A->p remains shallow on output if shallow on
15 // input). The A->x and A->i content is not changed; it remains in whatever
16 // shallow/non-shallow state that it had on input).
17
18 // If an out-of-memory condition occurs, all content of the matrix is cleared.
19
20 // If the input matrix A is sparse, bitmap or full, it is unchanged.
21
22 #include "GB.h"
23
24 GB_PUBLIC // accessed by the MATLAB tests in GraphBLAS/Test only
GB_convert_hyper_to_sparse(GrB_Matrix A,GB_Context Context)25 GrB_Info GB_convert_hyper_to_sparse // convert hypersparse to sparse
26 (
27 GrB_Matrix A, // matrix to convert to non-hypersparse
28 GB_Context Context
29 )
30 {
31
32 //--------------------------------------------------------------------------
33 // check inputs
34 //--------------------------------------------------------------------------
35
36 ASSERT_MATRIX_OK (A, "A being converted from hyper to sparse", GB0) ;
37 ASSERT (GB_ZOMBIES_OK (A)) ;
38 ASSERT (GB_JUMBLED_OK (A)) ;
39 ASSERT (GB_PENDING_OK (A)) ;
40
41 //--------------------------------------------------------------------------
42 // convert A from hypersparse to sparse
43 //--------------------------------------------------------------------------
44
45 if (GB_IS_HYPERSPARSE (A))
46 {
47
48 //----------------------------------------------------------------------
49 // determine the number of threads to use
50 //----------------------------------------------------------------------
51
52 GBURBLE ("(hyper to sparse) ") ;
53 int64_t n = A->vdim ;
54
55 GB_GET_NTHREADS_MAX (nthreads_max, chunk, Context) ;
56 int nthreads = GB_nthreads (n, chunk, nthreads_max) ;
57 int ntasks = (nthreads == 1) ? 1 : (8 * nthreads) ;
58 ntasks = GB_IMIN (ntasks, n) ;
59 ntasks = GB_IMAX (ntasks, 1) ;
60
61 //----------------------------------------------------------------------
62 // allocate the new Ap array, of size n+1
63 //----------------------------------------------------------------------
64
65 int64_t *restrict Ap_new = NULL ; size_t Ap_new_size = 0 ;
66 Ap_new = GB_MALLOC (n+1, int64_t, &Ap_new_size) ;
67 if (Ap_new == NULL)
68 {
69 // out of memory
70 return (GrB_OUT_OF_MEMORY) ;
71 }
72
73 #ifdef GB_DEBUG
74 // to ensure all values of Ap_new are assigned below.
75 for (int64_t j = 0 ; j <= n ; j++) Ap_new [j] = -99999 ;
76 #endif
77
78 //----------------------------------------------------------------------
79 // get the old hyperlist
80 //----------------------------------------------------------------------
81
82 int64_t nvec = A->nvec ; // # of vectors in Ah_old
83 int64_t *restrict Ap_old = A->p ; // size nvec+1
84 int64_t *restrict Ah_old = A->h ; // size nvec
85 int64_t nvec_nonempty = 0 ; // recompute A->nvec_nonempty
86 int64_t anz = GB_NNZ (A) ;
87
88 //----------------------------------------------------------------------
89 // construct the new vector pointers
90 //----------------------------------------------------------------------
91
92 int tid ;
93 #pragma omp parallel for num_threads(nthreads) schedule(dynamic,1) \
94 reduction(+:nvec_nonempty)
95 for (tid = 0 ; tid < ntasks ; tid++)
96 {
97 int64_t jstart, jend, my_nvec_nonempty = 0 ;
98 GB_PARTITION (jstart, jend, n, tid, ntasks) ;
99 ASSERT (0 <= jstart && jstart <= jend && jend <= n) ;
100
101 // task tid computes Ap_new [jstart:jend-1] from Ap_old, Ah_old.
102
103 // GB_SPLIT_BINARY_SEARCH of Ah_old [0..nvec-1] for jstart:
104 // If found is true then Ah_old [k] == jstart.
105 // If found is false, and nvec > 0 then
106 // Ah_old [0 ... k-1] < jstart < Ah_old [k ... nvec-1]
107 // Whether or not i is found, if nvec > 0
108 // Ah_old [0 ... k-1] < jstart <= Ah_old [k ... nvec-1]
109 // If nvec == 0, then k == 0 and found will be false. In this
110 // case, jstart cannot be compared with any content of Ah_old,
111 // since Ah_old is completely empty (Ah_old [0] is invalid).
112
113 int64_t k = 0, pright = nvec-1 ;
114 bool found ;
115 GB_SPLIT_BINARY_SEARCH (jstart, Ah_old, k, pright, found) ;
116 ASSERT (k >= 0 && k <= nvec) ;
117 ASSERT (GB_IMPLIES (nvec == 0, !found && k == 0)) ;
118 ASSERT (GB_IMPLIES (found, jstart == Ah_old [k])) ;
119 ASSERT (GB_IMPLIES (!found && k < nvec, jstart < Ah_old [k])) ;
120
121 // Let jk = Ah_old [k], jlast = Ah_old [k-1], and pk = Ah_old [k].
122 // Then Ap_new [jlast+1:jk] must be set to pk. This must be done
123 // for all k = 0:nvec-1. In addition, the last vector k=nvec-1
124 // must be terminated by setting Ap_new [jk+1:n-1] to Ap_old [nvec].
125 // A task owns the kth vector if jk is in jstart:jend-1, inclusive.
126 // It counts all non-empty vectors that it owns. However, the task
127 // must also set Ap_new [...] = pk for any jlast+1:jk that overlaps
128 // jstart:jend-1, even if it does not own that particular vector k.
129 // This happens only at the tail end of jstart:jend-1.
130
131 int64_t jlast = (k == 0) ? (-1) : Ah_old [k-1] ;
132 jlast = GB_IMAX (jstart-1, jlast) ;
133
134 bool done = false ;
135
136 for ( ; k <= nvec && !done ; k++)
137 {
138
139 //--------------------------------------------------------------
140 // get the kth vector in Ah_old, which is vector index jk.
141 //--------------------------------------------------------------
142
143 int64_t jk = (k < nvec) ? Ah_old [k] : n ;
144 int64_t pk = (k < nvec) ? Ap_old [k] : anz ;
145
146 //--------------------------------------------------------------
147 // determine if this task owns jk
148 //--------------------------------------------------------------
149
150 int64_t jfin ;
151 if (jk >= jend)
152 {
153 // This is the last iteration for this task. This task
154 // does not own the kth vector. However, it does own the
155 // vector indices jlast+1:jend-1, and these vectors must
156 // be handled by this task.
157 jfin = jend - 1 ;
158 done = true ;
159 }
160 else
161 {
162 // This task owns the kth vector, which is vector index jk.
163 // Ap must be set to pk for all vector indices jlast+1:jk.
164 jfin = jk ;
165 ASSERT (k >= 0 && k < nvec && nvec > 0) ;
166 if (pk < Ap_old [k+1]) my_nvec_nonempty++ ;
167 }
168
169 //--------------------------------------------------------------
170 // set Ap_new for this vector
171 //--------------------------------------------------------------
172
173 // Ap_new [jlast+1:jk] must be set to pk. This tasks handles
174 // the intersection of jlast+1:jk with jstart:jend-1.
175
176 for (int64_t j = jlast+1 ; j <= jfin ; j++)
177 {
178 Ap_new [j] = pk ;
179 }
180
181 //--------------------------------------------------------------
182 // keep track of the prior vector index
183 //--------------------------------------------------------------
184
185 jlast = jk ;
186 }
187 nvec_nonempty += my_nvec_nonempty ;
188
189 //------------------------------------------------------------------
190 // no task owns Ap_new [n] so it is set by the last task
191 //------------------------------------------------------------------
192
193 if (tid == ntasks-1)
194 {
195 ASSERT (jend == n) ;
196 Ap_new [n] = anz ;
197 }
198 }
199
200 // free the old A->p and A->h hyperlist content.
201 // this clears A->nvec_nonempty so it must be restored below.
202 GB_ph_free (A) ;
203
204 // transplant the new vector pointers; matrix is no longer hypersparse
205 A->p = Ap_new ; A->p_size = Ap_new_size ;
206 A->h = NULL ;
207 A->nvec = n ;
208 A->nvec_nonempty = nvec_nonempty ;
209 A->plen = n ;
210 A->p_shallow = false ;
211 A->h_shallow = false ;
212 A->magic = GB_MAGIC ;
213 ASSERT (anz == GB_NNZ (A)) ;
214
215 //----------------------------------------------------------------------
216 // A is now sparse
217 //----------------------------------------------------------------------
218
219 ASSERT (GB_IS_SPARSE (A)) ;
220 }
221
222 //--------------------------------------------------------------------------
223 // A is now in sparse form (or left as full or bitmap)
224 //--------------------------------------------------------------------------
225
226 ASSERT_MATRIX_OK (A, "A converted to sparse (or left as-is)", GB0) ;
227 ASSERT (!GB_IS_HYPERSPARSE (A)) ;
228 ASSERT (GB_ZOMBIES_OK (A)) ;
229 ASSERT (GB_JUMBLED_OK (A)) ;
230 ASSERT (GB_PENDING_OK (A)) ;
231 return (GrB_SUCCESS) ;
232 }
233
234