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
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
|
//------------------------------------------------------------------------------
// GB_convert_bitmap_worker: construct triplets or CSC/CSR from bitmap
//------------------------------------------------------------------------------
// SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2022, All Rights Reserved.
// SPDX-License-Identifier: Apache-2.0
//------------------------------------------------------------------------------
// If A is iso and Ax_new is not NULL, the iso scalar is expanded into the
// non-iso array Ax_new. Otherwise, if Ax_new and Ax are NULL then no values
// are extracted.
// TODO allow this function to do typecasting. Create 169 different versions
// for all 13x13 versions. Use this as part of Method 24, C=A assignment.
// Can also use typecasting for GB_Matrix_diag.
#include "GB.h"
#include "GB_partition.h"
#include "GB_unused.h"
GrB_Info GB_convert_bitmap_worker // extract CSC/CSR or triplets from bitmap
(
// outputs:
int64_t *restrict Ap, // vector pointers for CSC/CSR form
int64_t *restrict Ai, // indices for CSC/CSR or triplet form
int64_t *restrict Aj, // vector indices for triplet form
GB_void *restrict Ax_new, // values for CSC/CSR or triplet form
int64_t *anvec_nonempty, // # of non-empty vectors
// inputs: not modified
const GrB_Matrix A, // matrix to extract; not modified
GB_Context Context
)
{
//--------------------------------------------------------------------------
// check inputs
//--------------------------------------------------------------------------
ASSERT (GB_IS_BITMAP (A)) ;
ASSERT (Ap != NULL) ; // must be provided on input, size avdim+1
int64_t *restrict W = NULL ; size_t W_size = 0 ;
const int64_t avdim = A->vdim ;
const int64_t avlen = A->vlen ;
const size_t asize = A->type->size ;
//--------------------------------------------------------------------------
// count the entries in each vector
//--------------------------------------------------------------------------
const int8_t *restrict Ab = A->b ;
GB_GET_NTHREADS_MAX (nthreads_max, chunk, Context) ;
int nthreads = GB_nthreads (avlen*avdim, chunk, nthreads_max) ;
bool by_vector = (nthreads <= avdim) ;
if (by_vector)
{
//----------------------------------------------------------------------
// compute all vectors in parallel (no workspace)
//----------------------------------------------------------------------
int64_t j ;
#pragma omp parallel for num_threads(nthreads) schedule(static)
for (j = 0 ; j < avdim ; j++)
{
// ajnz = nnz (A (:,j))
int64_t ajnz = 0 ;
int64_t pA_start = j * avlen ;
for (int64_t i = 0 ; i < avlen ; i++)
{
// see if A(i,j) is present in the bitmap
int64_t p = i + pA_start ;
ajnz += Ab [p] ;
ASSERT (Ab [p] == 0 || Ab [p] == 1) ;
}
Ap [j] = ajnz ;
}
}
else
{
//----------------------------------------------------------------------
// compute blocks of rows in parallel
//----------------------------------------------------------------------
// allocate one row of W per thread, each row of length avdim
W = GB_MALLOC_WORK (nthreads * avdim, int64_t, &W_size) ;
if (W == NULL)
{
// out of memory
return (GrB_OUT_OF_MEMORY) ;
}
int taskid ;
#pragma omp parallel for num_threads(nthreads) schedule(static)
for (taskid = 0 ; taskid < nthreads ; taskid++)
{
int64_t *restrict Wtask = W + taskid * avdim ;
int64_t istart, iend ;
GB_PARTITION (istart, iend, avlen, taskid, nthreads) ;
for (int64_t j = 0 ; j < avdim ; j++)
{
// ajnz = nnz (A (istart:iend-1,j))
int64_t ajnz = 0 ;
int64_t pA_start = j * avlen ;
for (int64_t i = istart ; i < iend ; i++)
{
// see if A(i,j) is present in the bitmap
int64_t p = i + pA_start ;
ajnz += Ab [p] ;
ASSERT (Ab [p] == 0 || Ab [p] == 1) ;
}
Wtask [j] = ajnz ;
}
}
// cumulative sum to compute nnz(A(:,j)) for each vector j
int64_t j ;
#pragma omp parallel for num_threads(nthreads) schedule(static)
for (j = 0 ; j < avdim ; j++)
{
int64_t ajnz = 0 ;
for (int taskid = 0 ; taskid < nthreads ; taskid++)
{
int64_t *restrict Wtask = W + taskid * avdim ;
int64_t c = Wtask [j] ;
Wtask [j] = ajnz ;
ajnz += c ;
}
Ap [j] = ajnz ;
}
}
//--------------------------------------------------------------------------
// cumulative sum of Ap
//--------------------------------------------------------------------------
int nth = GB_nthreads (avdim, chunk, nthreads_max) ;
GB_cumsum (Ap, avdim, anvec_nonempty, nth, Context) ;
ASSERT (Ap [avdim] == A->nvals) ;
//--------------------------------------------------------------------------
// gather the pattern and values from the bitmap
//--------------------------------------------------------------------------
// TODO: add type-specific versions for built-in types
const GB_void *restrict Ax = (GB_void *) (A->x) ;
const bool A_iso = A->iso ;
const bool numeric = (Ax_new != NULL && Ax != NULL) ;
if (by_vector)
{
//----------------------------------------------------------------------
// construct all vectors in parallel (no workspace)
//----------------------------------------------------------------------
int64_t j ;
#pragma omp parallel for num_threads(nthreads) schedule(static)
for (j = 0 ; j < avdim ; j++)
{
// gather from the bitmap into the new A (:,j)
int64_t pnew = Ap [j] ;
int64_t pA_start = j * avlen ;
for (int64_t i = 0 ; i < avlen ; i++)
{
int64_t p = i + pA_start ;
if (Ab [p])
{
// A(i,j) is in the bitmap
if (Ai != NULL) Ai [pnew] = i ;
if (Aj != NULL) Aj [pnew] = j ;
if (numeric)
{
// Ax_new [pnew] = Ax [p])
memcpy (Ax_new +(pnew)*asize,
Ax +(A_iso ? 0:(p)*asize), asize) ;
}
pnew++ ;
}
}
ASSERT (pnew == Ap [j+1]) ;
}
}
else
{
//----------------------------------------------------------------------
// compute blocks of rows in parallel
//----------------------------------------------------------------------
int taskid ;
#pragma omp parallel for num_threads(nthreads) schedule(static)
for (taskid = 0 ; taskid < nthreads ; taskid++)
{
int64_t *restrict Wtask = W + taskid * avdim ;
int64_t istart, iend ;
GB_PARTITION (istart, iend, avlen, taskid, nthreads) ;
for (int64_t j = 0 ; j < avdim ; j++)
{
// gather from the bitmap into the new A (:,j)
int64_t pnew = Ap [j] + Wtask [j] ;
int64_t pA_start = j * avlen ;
for (int64_t i = istart ; i < iend ; i++)
{
// see if A(i,j) is present in the bitmap
int64_t p = i + pA_start ;
if (Ab [p])
{
// A(i,j) is in the bitmap
if (Ai != NULL) Ai [pnew] = i ;
if (Aj != NULL) Aj [pnew] = j ;
if (numeric)
{
// Ax_new [pnew] = Ax [p] ;
memcpy (Ax_new +(pnew)*asize,
Ax +(A_iso ? 0:(p)*asize), asize) ;
}
pnew++ ;
}
}
}
}
}
//--------------------------------------------------------------------------
// free workspace return result
//--------------------------------------------------------------------------
GB_FREE_WORK (&W, W_size) ;
return (GrB_SUCCESS) ;
}
|
