tatami
C++ API for different matrix representations
Loading...
Searching...
No Matches
convert_to_compressed_sparse.hpp
Go to the documentation of this file.
1#ifndef TATAMI_CONVERT_TO_COMPRESSED_SPARSE_H
2#define TATAMI_CONVERT_TO_COMPRESSED_SPARSE_H
3
8
9#include <memory>
10#include <vector>
11
18namespace tatami {
19
23namespace convert_to_compressed_sparse_internal {
24
25template<typename Value_, typename Index_, typename Count_>
26void count_compressed_sparse_non_zeros_consistent(const tatami::Matrix<Value_, Index_>& matrix, Index_ primary, Index_ secondary, bool row, Count_* output, int threads) {
27 if (matrix.is_sparse()) {
28 Options opt;
29 opt.sparse_extract_value = false;
30 opt.sparse_extract_index = false;
31 opt.sparse_ordered_index = false;
32
33 parallelize([&](int, Index_ start, Index_ length) -> void {
34 auto wrk = consecutive_extractor<true>(matrix, row, start, length, opt);
35 for (Index_ x = 0; x < length; ++x) {
36 auto range = wrk->fetch(NULL, NULL);
37 output[start + x] = range.number;
38 }
39 }, primary, threads);
40
41 } else {
42 parallelize([&](int, Index_ start, Index_ length) -> void {
43 std::vector<Value_> buffer_v(secondary);
44 auto wrk = consecutive_extractor<false>(matrix, row, start, length);
45 for (Index_ p = start, pe = start + length; p < pe; ++p) {
46 auto ptr = wrk->fetch(buffer_v.data());
47 Count_ count = 0;
48 for (Index_ s = 0; s < secondary; ++s) {
49 count += (ptr[s] != 0);
50 }
51 output[p] = count;
52 }
53 }, primary, threads);
54 }
55}
56
57template<typename Value_, typename Index_, typename Count_>
58void count_compressed_sparse_non_zeros_inconsistent(const tatami::Matrix<Value_, Index_>& matrix, Index_ primary, Index_ secondary, bool row, Count_* output, int threads) {
59 std::vector<std::vector<Count_> > nz_counts(threads - 1);
60 for (auto& x : nz_counts) {
61 x.resize(primary);
62 }
63
64 if (matrix.is_sparse()) {
65 Options opt;
66 opt.sparse_extract_value = false;
67 opt.sparse_ordered_index = false;
68
69 parallelize([&](int t, Index_ start, Index_ length) -> void {
70 std::vector<Index_> buffer_i(primary);
71 auto wrk = consecutive_extractor<true>(matrix, !row, start, length, opt);
72 auto my_counts = (t > 0 ? nz_counts[t - 1].data() : output);
73
74 for (Index_ x = 0; x < length; ++x) {
75 auto range = wrk->fetch(NULL, buffer_i.data());
76 for (Index_ i = 0; i < range.number; ++i) {
77 ++my_counts[range.index[i]];
78 }
79 }
80 }, secondary, threads);
81
82 } else {
83 parallelize([&](int t, Index_ start, Index_ length) -> void {
84 auto wrk = consecutive_extractor<false>(matrix, !row, start, length);
85 std::vector<Value_> buffer_v(primary);
86 auto my_counts = (t > 0 ? nz_counts[t - 1].data() : output);
87
88 for (Index_ x = 0; x < length; ++x) {
89 auto ptr = wrk->fetch(buffer_v.data());
90 for (Index_ p = 0; p < primary; ++p) {
91 my_counts[p] += (ptr[p] != 0);
92 }
93 }
94 }, secondary, threads);
95 }
96
97 for (auto& y : nz_counts) {
98 for (Index_ p = 0; p < primary; ++p) {
99 output[p] += y[p];
100 }
101 }
102}
103
104template<typename InputValue_, typename InputIndex_, typename Pointer_, typename StoredValue_, typename StoredIndex_>
105void fill_compressed_sparse_matrix_consistent(
107 InputIndex_ primary,
108 InputIndex_ secondary,
109 bool row,
110 const Pointer_* pointers,
111 StoredValue_* output_value,
112 StoredIndex_* output_index,
113 int threads)
114{
115 if (matrix.is_sparse()) {
116 Options opt;
117 opt.sparse_ordered_index = false;
118
119 parallelize([&](int, InputIndex_ start, InputIndex_ length) -> void {
120 std::vector<InputValue_> buffer_v(secondary);
121 std::vector<InputIndex_> buffer_i(secondary);
122 auto wrk = consecutive_extractor<true>(matrix, row, start, length, opt);
123
124 for (InputIndex_ p = start, pe = start + length; p < pe; ++p) {
125 // Resist the urge to `fetch()` straight into 'output_v'
126 // and 'output_i', as implementations may assume that they
127 // have the entire 'length' length to play with, and the
128 // output vectors only have whatever is allocated from the
129 // first pass (which might be nothing for an all-zero matrix).
130 auto range = wrk->fetch(buffer_v.data(), buffer_i.data());
131 auto offset = pointers[p];
132 std::copy_n(range.value, range.number, output_value + offset);
133 std::copy_n(range.index, range.number, output_index + offset);
134 }
135 }, primary, threads);
136
137 } else {
138 parallelize([&](int, InputIndex_ start, InputIndex_ length) -> void {
139 std::vector<InputValue_> buffer_v(secondary);
140 auto wrk = consecutive_extractor<false>(matrix, row, start, length);
141
142 for (InputIndex_ p = start, pe = start + length; p < pe; ++p) {
143 auto ptr = wrk->fetch(buffer_v.data());
144 auto offset = pointers[p];
145 for (InputIndex_ s = 0; s < secondary; ++s) {
146 auto val = ptr[s];
147 if (val != 0) {
148 output_value[offset] = val;
149 output_index[offset] = s;
150 ++offset;
151 }
152 }
153 }
154 }, primary, threads);
155 }
156}
157
158template<typename InputValue_, typename InputIndex_, typename Pointer_, typename StoredValue_, typename StoredIndex_>
159void fill_compressed_sparse_matrix_inconsistent(
161 InputIndex_ primary,
162 InputIndex_ secondary,
163 bool row,
164 const Pointer_* pointers,
165 StoredValue_* output_value,
166 StoredIndex_* output_index,
167 int threads)
168{
169 if (matrix.is_sparse()) {
170 Options opt;
171 opt.sparse_ordered_index = false;
172
173 parallelize([&](int, InputIndex_ start, InputIndex_ length) -> void {
174 std::vector<InputValue_> buffer_v(length);
175 std::vector<InputIndex_> buffer_i(length);
176 auto wrk = consecutive_extractor<true>(matrix, !row, static_cast<InputIndex_>(0), secondary, start, length, opt);
177 std::vector<Pointer_> offset_copy(pointers + start, pointers + start + length);
178
179 for (InputIndex_ x = 0; x < secondary; ++x) {
180 auto range = wrk->fetch(buffer_v.data(), buffer_i.data());
181 for (InputIndex_ i = 0; i < range.number; ++i) {
182 auto& pos = offset_copy[range.index[i] - start];
183 output_value[pos] = range.value[i];
184 output_index[pos] = x;
185 ++pos;
186 }
187 }
188 }, primary, threads);
189
190 } else {
191 parallelize([&](int, InputIndex_ start, InputIndex_ length) -> void {
192 std::vector<InputValue_> buffer_v(length);
193 auto wrk = consecutive_extractor<false>(matrix, !row, static_cast<InputIndex_>(0), secondary, start, length);
194 std::vector<Pointer_> offset_copy(pointers + start, pointers + start + length);
195
196 for (InputIndex_ x = 0; x < secondary; ++x) {
197 auto ptr = wrk->fetch(buffer_v.data());
198 for (InputIndex_ p = 0; p < length; ++p) {
199 auto val = ptr[p];
200 if (val != 0) {
201 auto& pos = offset_copy[p];
202 output_value[pos] = val;
203 output_index[pos] = x;
204 ++pos;
205 }
206 }
207 }
208 }, primary, threads);
209 }
210}
211
212}
226
242template<typename Value_, typename Index_, typename Count_>
244 Index_ NR = matrix.nrow();
245 Index_ NC = matrix.ncol();
246 Index_ primary = (row ? NR : NC);
247 Index_ secondary = (row ? NC : NR);
248 std::fill_n(output, primary, 0);
249
250 if (row == matrix.prefer_rows()) {
251 convert_to_compressed_sparse_internal::count_compressed_sparse_non_zeros_consistent(matrix, primary, secondary, row, output, options.num_threads);
252 } else {
253 convert_to_compressed_sparse_internal::count_compressed_sparse_non_zeros_inconsistent(matrix, primary, secondary, row, output, options.num_threads);
254 }
255}
256
266
286template<typename InputValue_, typename InputIndex_, typename Pointer_, typename StoredValue_, typename StoredIndex_>
289 bool row,
290 const Pointer_* pointers,
291 StoredValue_* output_value,
292 StoredIndex_* output_index,
294{
295 InputIndex_ NR = matrix.nrow();
296 InputIndex_ NC = matrix.ncol();
297 InputIndex_ primary = (row ? NR : NC);
298 InputIndex_ secondary = (row ? NC : NR);
299
300 if (row == matrix.prefer_rows()) {
301 convert_to_compressed_sparse_internal::fill_compressed_sparse_matrix_consistent(matrix, primary, secondary, row, pointers, output_value, output_index, options.num_threads);
302 } else {
303 convert_to_compressed_sparse_internal::fill_compressed_sparse_matrix_inconsistent(matrix, primary, secondary, row, pointers, output_value, output_index, options.num_threads);
304 }
305}
306
317template<typename Value_, typename Index_, typename Pointer_>
322 std::vector<Value_> value;
323
327 std::vector<Index_> index;
328
332 std::vector<Pointer_> pointers;
333};
334
350
368template<typename StoredValue_, typename StoredIndex_, typename StoredPointer_ = size_t, typename InputValue_, typename InputIndex_>
371 bool row,
373{
375 auto& output_v = output.value;
376 auto& output_i = output.index;
377 auto& output_p = output.pointers;
378
379 InputIndex_ NR = matrix.nrow();
380 InputIndex_ NC = matrix.ncol();
381 InputIndex_ primary = (row ? NR : NC);
382 InputIndex_ secondary = (row ? NC : NR);
383
384 if (!options.two_pass) {
385 // Doing a single fragmented run and then concatenating everything together.
387 matrix,
388 row,
389 [&]{
391 roptions.num_threads = options.num_threads;
392 return roptions;
393 }()
394 );
395 const auto& store_v = frag.value;
396 const auto& store_i = frag.index;
397
398 output_p.resize(static_cast<size_t>(primary) + 1);
399 for (InputIndex_ p = 0; p < primary; ++p) {
400 output_p[p + 1] = output_p[p] + store_v[p].size();
401 }
402
403 output_v.reserve(output_p.back());
404 output_i.reserve(output_p.back());
405 for (InputIndex_ p = 0; p < primary; ++p) {
406 output_v.insert(output_v.end(), store_v[p].begin(), store_v[p].end());
407 output_i.insert(output_i.end(), store_i[p].begin(), store_i[p].end());
408 }
409
410 } else if (row == matrix.prefer_rows()) {
411 // First pass to figure out how many non-zeros there are.
412 output_p.resize(static_cast<size_t>(primary) + 1);
413 convert_to_compressed_sparse_internal::count_compressed_sparse_non_zeros_consistent(matrix, primary, secondary, row, output_p.data() + 1, options.num_threads);
414 for (InputIndex_ i = 1; i <= primary; ++i) {
415 output_p[i] += output_p[i - 1];
416 }
417
418 // Second pass to actually fill our vectors.
419 output_v.resize(output_p.back());
420 output_i.resize(output_p.back());
421 convert_to_compressed_sparse_internal::fill_compressed_sparse_matrix_consistent(
422 matrix,
423 primary,
424 secondary,
425 row,
426 output_p.data(),
427 output_v.data(),
428 output_i.data(),
429 options.num_threads
430 );
431
432 } else {
433 // First pass to figure out how many non-zeros there are.
434 output_p.resize(static_cast<size_t>(primary) + 1);
435 convert_to_compressed_sparse_internal::count_compressed_sparse_non_zeros_inconsistent(matrix, primary, secondary, row, output_p.data() + 1, options.num_threads);
436 for (InputIndex_ i = 1; i <= primary; ++i) {
437 output_p[i] += output_p[i - 1];
438 }
439
440 // Second pass to actually fill our vectors.
441 output_v.resize(output_p.back());
442 output_i.resize(output_p.back());
443 convert_to_compressed_sparse_internal::fill_compressed_sparse_matrix_inconsistent(
444 matrix,
445 primary,
446 secondary,
447 row,
448 output_p.data(),
449 output_v.data(),
450 output_i.data(),
451 options.num_threads
452 );
453 }
454
455 return output;
456}
457
466 bool two_pass = false;
467
471 int num_threads = 1;
472};
473
491template<
492 typename Value_,
493 typename Index_,
494 typename StoredValue_ = Value_,
495 typename StoredIndex_ = Index_,
496 typename StoredPointer_ = size_t,
497 typename InputValue_,
498 typename InputIndex_
499>
500std::shared_ptr<Matrix<Value_, Index_> > convert_to_compressed_sparse(const Matrix<InputValue_, InputIndex_>& matrix, bool row, const ConvertToCompressedSparseOptions& options) {
502 matrix,
503 row,
504 [&]{
506 ropt.two_pass = options.two_pass;
507 ropt.num_threads = options.num_threads;
508 return ropt;
509 }()
510 );
511 return std::shared_ptr<Matrix<Value_, Index_> >(
513 Value_,
514 Index_,
515 std::vector<StoredValue_>,
516 std::vector<StoredIndex_>,
517 std::vector<StoredPointer_>
518 >(
519 matrix.nrow(),
520 matrix.ncol(),
521 std::move(comp.value),
522 std::move(comp.index),
523 std::move(comp.pointers),
524 row,
525 []{
526 CompressedSparseMatrixOptions copt;
527 copt.check = false; // no need for checks, as we guarantee correctness.
528 return copt;
529 }()
530 )
531 );
532}
533
537// Backwards compatbility.
538template<typename Value_, typename Index_, typename Count_>
539void count_compressed_sparse_non_zeros(const tatami::Matrix<Value_, Index_>* matrix, bool row, Count_* output, int threads) {
541 *matrix,
542 row,
543 output,
544 [&]{
545 CountCompressedSparseNonZerosOptions copt;
546 copt.num_threads = threads;
547 return copt;
548 }()
549 );
550}
551
552template<typename InputValue_, typename InputIndex_, typename Pointer_, typename StoredValue_, typename StoredIndex_>
554 bool row,
555 const Pointer_* pointers,
556 StoredValue_* output_value,
557 StoredIndex_* output_index,
558 int threads)
559{
561 *matrix,
562 row,
563 pointers,
564 output_value,
565 output_index,
566 [&]{
567 FillCompressedSparseContentsOptions fopt;
568 fopt.num_threads = threads;
569 return fopt;
570 }()
571 );
572}
573
574template<typename StoredValue_, typename StoredIndex_, typename StoredPointer_ = size_t, typename InputValue_, typename InputIndex_>
575CompressedSparseContents<StoredValue_, StoredIndex_, StoredPointer_> retrieve_compressed_sparse_contents(const Matrix<InputValue_, InputIndex_>* matrix, bool row, bool two_pass, int threads = 1) {
577 *matrix,
578 row,
579 [&]{
580 RetrieveCompressedSparseContentsOptions opt;
581 opt.two_pass = two_pass;
582 opt.num_threads = threads;
583 return opt;
584 }()
585 );
586}
587
588template<typename Value_ = double, typename Index_ = int, typename StoredValue_ = Value_, typename StoredIndex_ = Index_, typename InputValue_, typename InputIndex_>
589std::shared_ptr<Matrix<Value_, Index_> > convert_to_compressed_sparse(const Matrix<InputValue_, InputIndex_>* matrix, bool row, bool two_pass = false, int threads = 1) {
591 *matrix,
592 row,
593 [&]{
594 ConvertToCompressedSparseOptions opt;
595 opt.two_pass = two_pass;
596 opt.num_threads = threads;
597 return opt;
598 }()
599 );
600}
601
602template <bool row_, typename Value_, typename Index_, typename InputValue_, typename InputIndex_>
603CompressedSparseContents<Value_, Index_, size_t> retrieve_compressed_sparse_contents(const Matrix<InputValue_, InputIndex_>* matrix, bool two_pass, int threads = 1) {
604 return retrieve_compressed_sparse_contents<Value_, Index_>(matrix, row_, two_pass, threads);
605}
606
607template <bool row_, typename Value_, typename Index_, typename StoredValue_ = Value_, typename StoredIndex_ = Index_, typename InputValue_, typename InputIndex_>
608std::shared_ptr<Matrix<Value_, Index_> > convert_to_compressed_sparse(const Matrix<InputValue_, InputIndex_>* matrix, bool two_pass = false, int threads = 1) {
610}
615}
616
617#endif
Compressed sparse matrix representation.
Compressed sparse matrix representation.
Definition CompressedSparseMatrix.hpp:505
Virtual class for a matrix.
Definition Matrix.hpp:59
virtual Index_ ncol() const =0
virtual Index_ nrow() const =0
virtual bool prefer_rows() const =0
virtual bool is_sparse() const =0
Templated construction of a new consecutive extractor.
Convert a matrix into a fragmented sparse format.
Flexible representations for matrix data.
Definition Extractor.hpp:15
CompressedSparseContents< StoredValue_, StoredIndex_, StoredPointer_ > retrieve_compressed_sparse_contents(const Matrix< InputValue_, InputIndex_ > &matrix, bool row, const RetrieveCompressedSparseContentsOptions &options)
Definition convert_to_compressed_sparse.hpp:369
std::shared_ptr< Matrix< Value_, Index_ > > convert_to_compressed_sparse(const Matrix< InputValue_, InputIndex_ > &matrix, bool row, const ConvertToCompressedSparseOptions &options)
Definition convert_to_compressed_sparse.hpp:500
void parallelize(Function_ fun, Index_ tasks, int threads)
Definition parallelize.hpp:42
FragmentedSparseContents< StoredValue_, StoredIndex_ > retrieve_fragmented_sparse_contents(const Matrix< InputValue_, InputIndex_ > &matrix, bool row, const RetrieveFragmentedSparseContentsOptions &options)
Definition convert_to_fragmented_sparse.hpp:76
void fill_compressed_sparse_contents(const tatami::Matrix< InputValue_, InputIndex_ > &matrix, bool row, const Pointer_ *pointers, StoredValue_ *output_value, StoredIndex_ *output_index, const FillCompressedSparseContentsOptions &options)
Definition convert_to_compressed_sparse.hpp:287
void count_compressed_sparse_non_zeros(const tatami::Matrix< Value_, Index_ > &matrix, bool row, Count_ *output, const CountCompressedSparseNonZerosOptions &options)
Definition convert_to_compressed_sparse.hpp:243
Parallelized iteration over a tatami::Matrix.
Compressed sparse contents.
Definition convert_to_compressed_sparse.hpp:318
std::vector< Index_ > index
Definition convert_to_compressed_sparse.hpp:327
std::vector< Value_ > value
Definition convert_to_compressed_sparse.hpp:322
std::vector< Pointer_ > pointers
Definition convert_to_compressed_sparse.hpp:332
Options for convert_to_compressed_sparse().
Definition convert_to_compressed_sparse.hpp:461
bool two_pass
Definition convert_to_compressed_sparse.hpp:466
int num_threads
Definition convert_to_compressed_sparse.hpp:471
Options for count_compressed_sparse_non_zeros().
Definition convert_to_compressed_sparse.hpp:220
int num_threads
Definition convert_to_compressed_sparse.hpp:224
Options for fill_compressed_sparse_contents().
Definition convert_to_compressed_sparse.hpp:260
int num_threads
Definition convert_to_compressed_sparse.hpp:264
Options for retrieve_compressed_sparse_contents().
Definition convert_to_compressed_sparse.hpp:338
int num_threads
Definition convert_to_compressed_sparse.hpp:348
bool two_pass
Definition convert_to_compressed_sparse.hpp:343
Options for retrieve_fragmented_sparse_contents().
Definition convert_to_fragmented_sparse.hpp:56
int num_threads
Definition convert_to_fragmented_sparse.hpp:60