MVE - Multi-View Environment mve-devel
Loading...
Searching...
No Matches
ba_sparse_matrix.h
Go to the documentation of this file.
1/*
2 * Copyright (C) 2015, Simon Fuhrmann, Fabian Langguth
3 * TU Darmstadt - Graphics, Capture and Massively Parallel Computing
4 * All rights reserved.
5 *
6 * This software may be modified and distributed under the terms
7 * of the BSD 3-Clause license. See the LICENSE.txt file for details.
8 */
9
10#ifndef SFM_SPARSE_MATRIX_HEADER
11#define SFM_SPARSE_MATRIX_HEADER
12
13#include <thread>
14#include <stdexcept>
15#include <vector>
16#include <algorithm>
17
18#include "sfm/ba_dense_vector.h"
19#include "sfm/defines.h"
20
23
27template <typename T>
29{
30public:
32 struct Triplet
33 {
34 Triplet (void) = default;
35 Triplet (std::size_t row, std::size_t col, T const& value);
36
37 std::size_t row;
38 std::size_t col;
40 };
41
43 typedef std::vector<Triplet> Triplets;
44
45public:
46 SparseMatrix (void);
47 SparseMatrix (std::size_t rows, std::size_t cols);
48 void allocate (std::size_t rows, std::size_t cols);
49 void reserve (std::size_t num_elements);
50 void set_from_triplets (Triplets const& triplets);
51 void mult_diagonal (T const& factor);
52 void cwise_invert (void);
53 void column_nonzeros (std::size_t col, DenseVector<T>* vector) const;
54
55 SparseMatrix transpose (void) const;
56 SparseMatrix subtract (SparseMatrix const& rhs) const;
57 SparseMatrix multiply (SparseMatrix const& rhs) const;
58 SparseMatrix sequential_multiply (SparseMatrix const& rhs) const;
59 SparseMatrix parallel_multiply (SparseMatrix const& rhs) const;
60 DenseVector<T> multiply (DenseVector<T> const& rhs) const;
61 SparseMatrix diagonal_matrix (void) const;
62
63 std::size_t num_non_zero (void) const;
64 std::size_t num_rows (void) const;
65 std::size_t num_cols (void) const;
66 T* begin (void);
67 T* end (void);
68
69 void debug (void) const;
70
71private:
72 std::size_t rows;
73 std::size_t cols;
74 std::vector<T> values;
75 std::vector<std::size_t> outer;
76 std::vector<std::size_t> inner;
77};
78
81
82/* ------------------------ Implementation ------------------------ */
83
84#include <iostream>
85
88
89template <typename T>
91 std::size_t col, T const& value)
92 : row(row), col(col), value(value)
93{
94}
95
96/* --------------------------------------------------------------- */
97
98template <typename T>
100 : rows(0)
101 , cols(0)
102{
103}
104
105template <typename T>
106SparseMatrix<T>::SparseMatrix (std::size_t rows, std::size_t cols)
107{
108 this->allocate(rows, cols);
109}
110
111template <typename T>
112void
113SparseMatrix<T>::allocate (std::size_t rows, std::size_t cols)
114{
115 this->rows = rows;
116 this->cols = cols;
117 this->values.clear();
118 this->outer.clear();
119 this->inner.clear();
120 this->outer.resize(cols + 1, 0);
121}
122
123template <typename T>
124void
125SparseMatrix<T>::reserve (std::size_t num_elements)
126{
127 this->inner.reserve(num_elements);
128 this->values.reserve(num_elements);
129}
130
131template <typename T>
132void
134{
135 /* Create a temporary transposed matrix */
136 SparseMatrix<T> transposed(this->cols, this->rows);
137 transposed.values.resize(triplets.size());
138 transposed.inner.resize(triplets.size());
139
140 /* Initialize outer indices with amount of inner values. */
141 for (std::size_t i = 0; i < triplets.size(); ++i)
142 transposed.outer[triplets[i].row]++;
143
144 /* Convert amounts to indices with prefix sum. */
145 std::size_t sum = 0;
146 std::vector<std::size_t> scratch(transposed.outer.size());
147 for (std::size_t i = 0; i < transposed.outer.size(); ++i)
148 {
149 std::size_t const temp = transposed.outer[i];
150 transposed.outer[i] = sum;
151 scratch[i] = sum;
152 sum += temp;
153 }
154
155 /* Add triplets, inner indices are unsorted. */
156 for (std::size_t i = 0; i < triplets.size(); ++i)
157 {
158 Triplet const& t = triplets[i];
159 std::size_t pos = scratch[t.row]++;
160 transposed.values[pos] = t.value;
161 transposed.inner[pos] = t.col;
162 }
163
164 /* Transpose matrix, implicit sorting of inner indices. */
165 *this = transposed.transpose();
166}
167
168template <typename T>
171{
172 SparseMatrix ret(this->cols, this->rows);
173 ret.values.resize(this->num_non_zero());
174 ret.inner.resize(this->num_non_zero());
175
176 /* Compute inner sizes of transposed matrix. */
177 for(std::size_t i = 0; i < this->inner.size(); ++i)
178 ret.outer[this->inner[i]] += 1;
179
180 /* Compute outer sizes of transposed matrix with prefix sum. */
181 std::size_t sum = 0;
182 std::vector<std::size_t> scratch(ret.outer.size());
183 for (std::size_t i = 0; i < ret.outer.size(); ++i)
184 {
185 std::size_t const temp = ret.outer[i];
186 ret.outer[i] = sum;
187 scratch[i] = sum;
188 sum += temp;
189 }
190
191 /* Write inner indices and values of transposed matrix. */
192 for (std::size_t i = 0; i < this->outer.size() - 1; ++i)
193 for (std::size_t j = this->outer[i]; j < this->outer[i + 1]; ++j)
194 {
195 std::size_t pos = scratch[this->inner[j]]++;
196 ret.inner[pos] = i;
197 ret.values[pos] = this->values[j];
198 }
199
200 return ret;
201}
202
203template <typename T>
206{
207 if (this->rows != rhs.rows || this->cols != rhs.cols)
208 throw std::invalid_argument("Incompatible matrix dimensions");
209
210 SparseMatrix ret(this->rows, this->cols);
211 ret.reserve(this->num_non_zero() + rhs.num_non_zero());
212
213 std::size_t num_outer = this->outer.size() - 1;
214 for (std::size_t outer = 0; outer < num_outer; ++outer)
215 {
216 ret.outer[outer] = ret.values.size();
217
218 std::size_t i1 = this->outer[outer];
219 std::size_t i2 = rhs.outer[outer];
220 std::size_t const i1_end = this->outer[outer + 1];
221 std::size_t const i2_end = rhs.outer[outer + 1];
222 while (i1 < i1_end || i2 < i2_end)
223 {
224 if (i1 >= i1_end)
225 {
226 ret.values.push_back(-rhs.values[i2]);
227 ret.inner.push_back(rhs.inner[i2]);
228 i2 += 1;
229 continue;
230 }
231 if (i2 >= i2_end)
232 {
233 ret.values.push_back(this->values[i1]);
234 ret.inner.push_back(this->inner[i1]);
235 i1 += 1;
236 continue;
237 }
238
239 std::size_t id1 = this->inner[i1];
240 std::size_t id2 = rhs.inner[i2];
241
242 if (id1 < id2)
243 ret.values.push_back(this->values[i1]);
244 else if (id2 < id1)
245 ret.values.push_back(-rhs.values[i2]);
246 else
247 ret.values.push_back(this->values[i1] - rhs.values[i2]);
248
249 i1 += static_cast<std::size_t>(id1 <= id2);
250 i2 += static_cast<std::size_t>(id2 <= id1);
251 ret.inner.push_back(std::min(id1, id2));
252 }
253 }
254 ret.outer.back() = ret.values.size();
255
256 return ret;
257}
258
259template <typename T>
262{
263#ifdef _OPENMP
264 return this->parallel_multiply(rhs);
265#else
266 return this->sequential_multiply(rhs);
267#endif
268}
269
270template <typename T>
273{
274 if (this->cols != rhs.rows)
275 throw std::invalid_argument("Incompatible matrix dimensions");
276
277 SparseMatrix ret(this->rows, rhs.cols);
278 ret.reserve(this->num_non_zero() + rhs.num_non_zero());
279
280 /* Matrix-matrix multiplication. */
281 std::vector<T> ret_col(ret.rows, T(0));
282 std::vector<bool> ret_nonzero(ret.rows, false);
283 for (std::size_t col = 0; col < ret.cols; ++col)
284 {
285 ret.outer[col] = ret.values.size();
286
287 std::fill(ret_col.begin(), ret_col.end(), T(0));
288 std::fill(ret_nonzero.begin(), ret_nonzero.end(), false);
289 std::size_t rhs_col_begin = rhs.outer[col];
290 std::size_t rhs_col_end = rhs.outer[col + 1];
291 for (std::size_t i = rhs_col_begin; i < rhs_col_end; ++i)
292 {
293 T const& rhs_col_value = rhs.values[i];
294 std::size_t const lhs_col = rhs.inner[i];
295 std::size_t const lhs_col_begin = this->outer[lhs_col];
296 std::size_t const lhs_col_end = this->outer[lhs_col + 1];
297
298 for (std::size_t j = lhs_col_begin; j < lhs_col_end; ++j)
299 {
300 std::size_t const id = this->inner[j];
301 ret_col[id] += this->values[j] * rhs_col_value;
302 ret_nonzero[id] = true;
303 }
304 }
305 for (std::size_t i = 0; i < ret.rows; ++i)
306 if (ret_nonzero[i])
307 {
308 ret.inner.push_back(i);
309 ret.values.push_back(ret_col[i]);
310 }
311 }
312 ret.outer[ret.cols] = ret.values.size();
313
314 return ret;
315}
316
317template <typename T>
320{
321 if (this->cols != rhs.rows)
322 throw std::invalid_argument("Incompatible matrix dimensions");
323
324 std::size_t nnz = this->num_non_zero() + rhs.num_non_zero();
325 SparseMatrix ret(this->rows, rhs.cols);
326 ret.reserve(nnz);
327 std::fill(ret.outer.begin(), ret.outer.end(), 0);
328
329 std::size_t const chunk_size = 64;
330 std::size_t const num_chunks = ret.cols / chunk_size
331 + (ret.cols % chunk_size != 0);
332 std::size_t const max_threads = std::max(1u,
333 std::thread::hardware_concurrency());
334 std::size_t const num_threads = std::min(num_chunks, max_threads);
335
336#pragma omp parallel num_threads(num_threads)
337 {
338 /* Matrix-matrix multiplication. */
339 std::vector<T> ret_col(ret.rows, T(0));
340 std::vector<bool> ret_nonzero(ret.rows, false);
341 std::vector<T> thread_values;
342 thread_values.reserve(nnz / num_chunks);
343 std::vector<std::size_t> thread_inner;
344 thread_inner.reserve(nnz / num_chunks);
345
346#pragma omp for ordered schedule(static, 1)
347 for (std::size_t chunk = 0; chunk < num_chunks; ++chunk)
348 {
349 thread_inner.clear();
350 thread_values.clear();
351
352 std::size_t const begin = chunk * chunk_size;
353 std::size_t const end = std::min(begin + chunk_size, ret.cols);
354 for (std::size_t col = begin; col < end; ++col)
355 {
356 std::fill(ret_col.begin(), ret_col.end(), T(0));
357 std::fill(ret_nonzero.begin(), ret_nonzero.end(), false);
358 std::size_t const rhs_col_begin = rhs.outer[col];
359 std::size_t const rhs_col_end = rhs.outer[col + 1];
360 for (std::size_t i = rhs_col_begin; i < rhs_col_end; ++i)
361 {
362 T const& rhs_col_value = rhs.values[i];
363 std::size_t const lhs_col = rhs.inner[i];
364 std::size_t const lhs_col_begin = this->outer[lhs_col];
365 std::size_t const lhs_col_end = this->outer[lhs_col + 1];
366
367 for (std::size_t j = lhs_col_begin; j < lhs_col_end; ++j)
368 {
369 std::size_t const id = this->inner[j];
370 ret_col[id] += this->values[j] * rhs_col_value;
371 ret_nonzero[id] = true;
372 }
373 }
374 for (std::size_t i = 0; i < ret.rows; ++i)
375 if (ret_nonzero[i])
376 {
377 ret.outer[col + 1] += 1;
378 thread_inner.push_back(i);
379 thread_values.push_back(ret_col[i]);
380 }
381 }
382
383#pragma omp ordered
384 {
385 ret.inner.insert(ret.inner.end(),
386 thread_inner.begin(), thread_inner.end());
387 ret.values.insert(ret.values.end(),
388 thread_values.begin(), thread_values.end());
389 }
390 }
391 }
392
393 for (std::size_t col = 0; col < ret.cols; ++col)
394 ret.outer[col + 1] += ret.outer[col];
395
396 return ret;
397}
398
399template<typename T>
402{
403 if (rhs.size() != this->cols)
404 throw std::invalid_argument("Incompatible dimensions");
405
406 DenseVector<T> ret(this->rows, T(0));
407 for (std::size_t i = 0; i < this->cols; ++i)
408 for (std::size_t id = this->outer[i]; id < this->outer[i + 1]; ++id)
409 ret[this->inner[id]] += this->values[id] * rhs[i];
410 return ret;
411}
412
413template<typename T>
416{
417 std::size_t const diag_size = std::min(this->rows, this->cols);
418 SparseMatrix ret(diag_size, diag_size);
419 ret.reserve(diag_size);
420 for (std::size_t i = 0; i < diag_size; ++i)
421 {
422 ret.outer[i] = ret.values.size();
423 for (std::size_t j = this->outer[i]; j < this->outer[i + 1]; ++j)
424 if (this->inner[j] == i)
425 {
426 ret.inner.push_back(i);
427 ret.values.push_back(this->values[j]);
428 }
429 else if (this->inner[j] > i)
430 break;
431 }
432 ret.outer[diag_size] = ret.values.size();
433 return ret;
434}
435
436template<typename T>
437void
439{
440 for (std::size_t i = 0; i < this->outer.size() - 1; ++i)
441 for (std::size_t j = this->outer[i]; j < this->outer[i + 1]; ++j)
442 {
443 if (this->inner[j] == i)
444 this->values[j] *= factor;
445 if (this->inner[j] >= i)
446 break;
447 }
448}
449
450template<typename T>
451void
453{
454 for (std::size_t i = 0; i < this->values.size(); ++i)
455 this->values[i] = T(1) / this->values[i];
456}
457
458template<typename T>
459void
460SparseMatrix<T>::column_nonzeros (std::size_t col, DenseVector<T>* vector) const
461{
462 std::size_t const start = this->outer[col];
463 std::size_t const end = this->outer[col + 1];
464 vector->resize(end - start);
465 for (std::size_t row = start, i = 0; row < end; ++row, ++i)
466 vector->at(i) = this->values[row];
467}
468
469template<typename T>
470inline std::size_t
472{
473 return this->values.size();
474}
475
476template<typename T>
477inline std::size_t
479{
480 return this->rows;
481}
482
483template<typename T>
484inline std::size_t
486{
487 return this->cols;
488}
489
490template<typename T>
491inline T*
493{
494 return this->values.data();
495}
496
497template<typename T>
498inline T*
500{
501 return this->values.data() + this->values.size();
502}
503
504template<typename T>
505void
507{
508 std::cout << "SparseMatrix ("
509 << this->rows << " rows, " << this->cols << " cols, "
510 << this->num_non_zero() << " values)" << std::endl;
511 std::cout << " Values:";
512 for (std::size_t i = 0; i < this->values.size(); ++i)
513 std::cout << " " << this->values[i];
514 std::cout << std::endl << " Inner:";
515 for (std::size_t i = 0; i < this->inner.size(); ++i)
516 std::cout << " " << this->inner[i];
517 std::cout << std::endl << " Outer:";
518 for (std::size_t i = 0; i < this->outer.size(); ++i)
519 std::cout << " " << this->outer[i];
520 std::cout << std::endl;
521}
522
525
526#endif // SFM_SPARSE_MATRIX_HEADER
527
T & at(std::size_t index)
void resize(std::size_t size, T const &value=T(0))
std::size_t size(void) const
Sparse matrix class in Yale format for column-major matrices.
void column_nonzeros(std::size_t col, DenseVector< T > *vector) const
SparseMatrix diagonal_matrix(void) const
void mult_diagonal(T const &factor)
std::size_t num_rows(void) const
SparseMatrix parallel_multiply(SparseMatrix const &rhs) const
SparseMatrix transpose(void) const
void set_from_triplets(Triplets const &triplets)
std::vector< Triplet > Triplets
List of triplets.
SparseMatrix sequential_multiply(SparseMatrix const &rhs) const
SparseMatrix subtract(SparseMatrix const &rhs) const
void reserve(std::size_t num_elements)
std::size_t num_cols(void) const
SparseMatrix multiply(SparseMatrix const &rhs) const
std::size_t num_non_zero(void) const
void allocate(std::size_t rows, std::size_t cols)
#define SFM_BA_NAMESPACE_BEGIN
Definition defines.h:22
#define SFM_NAMESPACE_END
Definition defines.h:14
#define SFM_NAMESPACE_BEGIN
Definition defines.h:13
#define SFM_BA_NAMESPACE_END
Definition defines.h:23
Triplet with row/col index, and the actual value.