28 #ifndef SHARK_LINALG_BLAS_MATRIX_SPARSE_HPP 29 #define SHARK_LINALG_BLAS_MATRIX_SPARSE_HPP 36 template<
class T,
class I=std::
size_t>
40 typedef std::size_t size_type;
41 typedef std::ptrdiff_t difference_type;
43 typedef value_type scalar_type;
54 const_reference value()
const {
55 return const_cast<self_type const&
>(m_matrix)(m_i,m_j);
57 value_type& ref()
const {
59 index_type
const *start = m_matrix.
inner_indices() + m_matrix.m_rowStart[m_i];
60 index_type
const *end = m_matrix.
inner_indices() + m_matrix.m_rowEnd[m_i];
62 index_type
const *pos = std::lower_bound(start,end,m_j);
64 if (pos != end && *pos == m_j)
65 return m_matrix.m_values[(pos-start)+m_matrix.m_rowStart[m_i]];
71 pos-start + m_matrix.m_rowStart[m_i]
74 return *m_matrix.
set_element(posIter, m_j, m_matrix.m_zero);
81 m_matrix(m), m_i(i), m_j(j) {}
88 return ref() = other.value();
121 : m_size1(0), m_size2(0), m_nnz(0)
122 , m_rowStart(1,0), m_indices(0), m_values(0), m_zero(0) {}
125 : m_size1(size1), m_size2(size2), m_nnz(0)
126 , m_rowStart(size1 + 1,0)
128 , m_indices(non_zeros), m_values(non_zeros), m_zero(0) {}
132 : m_size1(e().
size1()), m_size2(e().
size2()), m_nnz(0)
133 , m_rowStart(e().
size1() + 1, 0)
134 , m_rowEnd(e().
size1(), 0)
135 , m_indices(non_zeros), m_values(non_zeros), m_zero(0) {
148 return m_values.size();
152 return m_rowStart[row+1] - m_rowStart[
row];
158 return m_rowEnd[
row] - m_rowStart[
row];
162 if(m_rowStart.empty())
164 return &m_rowStart[0];
174 return &m_indices[0];
183 if(m_rowStart.empty())
185 return &m_rowStart[0];
195 return &m_indices[0];
211 m_rowEnd[i] = m_rowStart[i]+non_zeros;
214 m_rowStart[
size1()] = m_rowEnd[i];
222 m_rowStart.resize(m_size1 + 1);
223 m_rowEnd.resize(m_size1);
224 std::fill(m_rowStart.begin(),m_rowStart.end(),0);
225 std::fill(m_rowEnd.begin(),m_rowEnd.end(),0);
230 m_indices.resize(non_zeros);
231 m_values.resize(non_zeros);
236 non_zeros =
std::min(m_size2,non_zeros);
238 std::size_t spaceDifference = non_zeros-
row_capacity(row);
241 if (spaceDifference >
nnz_capacity()-m_rowStart.back()) {
245 for (std::size_t i =
size1()-1; i !=
row; --i) {
247 value_type *valueRowEnd = this->
values() + m_rowEnd[i];
250 std::copy_backward(values,valueRowEnd, valueRowEnd+spaceDifference);
251 std::copy_backward(indices,indicesEnd, indicesEnd+spaceDifference);
252 m_rowStart[i]+=spaceDifference;
253 m_rowEnd[i]+=spaceDifference;
255 m_rowStart.back() +=spaceDifference;
272 index_type
const *pos = std::lower_bound(start,end,j);
274 if (pos != end && *pos == j)
275 return m_values[(pos-start)+m_rowStart[i]];
305 m_rowStart.swap(m.m_rowStart);
306 m_rowEnd.swap(m.m_rowEnd);
307 m_indices.swap(m.m_indices);
308 m_values.swap(m.m_values);
336 value_type* valuesi = a.
values() + a.m_rowStart[i];
337 value_type* valuesj = b.
values() + b.m_rowStart[j];
340 std::swap_ranges(indicesi,indicesi+nnzi,indicesj);
341 std::swap_ranges(valuesi, valuesi+nnzi,valuesj);
358 typedef compressed_storage_iterator<value_type, index_type const>
row_iterator;
366 const_row_iterator
row_end(std::size_t i)
const {
387 row_iterator
set_element(row_iterator pos, size_type index, value_type value) {
388 std::size_t
row = pos.row();
394 if (pos !=
row_end(row) && pos.index() == index) {
400 difference_type arrayPos = pos -
row_begin(row) + m_rowStart[
row];
408 m_values.begin()+arrayPos, m_values.begin() + m_rowEnd[
row],
409 m_values.begin() + m_rowEnd[
row] + 1
412 m_indices.begin()+arrayPos, m_indices.begin() + m_rowEnd[
row],
413 m_indices.begin() + m_rowEnd[
row] + 1
416 m_values[arrayPos] = value;
417 m_indices[arrayPos] = index;
427 std::size_t
row = start.row();
430 size_type rowEndPos = m_rowEnd[
row];
431 size_type rowStartPos = m_rowStart[
row];
432 difference_type rangeStartPos = start -
row_begin(row)+rowStartPos;
433 difference_type rangeEndPos = end -
row_begin(row)+rowStartPos;
434 difference_type rangeSize = end - start;
438 m_values.begin()+rangeEndPos,m_values.begin() + rowEndPos, m_values.begin() + rangeStartPos
441 m_indices.begin()+rangeEndPos,m_indices.begin() + rowEndPos , m_indices.begin() + rangeStartPos
443 m_rowEnd[
row] -= rangeSize;
451 row_iterator next = elem;
457 template<
class Archive>
459 ar &boost::serialization::make_nvp(
"outer_indices", m_rowStart);
460 ar &boost::serialization::make_nvp(
"outer_indices_end", m_rowEnd);
461 ar &boost::serialization::make_nvp(
"inner_indices", m_indices);
462 ar &boost::serialization::make_nvp(
"values", m_values);
469 std::vector<index_type> m_rowStart;
470 std::vector<index_type> m_rowEnd;
471 std::vector<index_type> m_indices;
472 std::vector<value_type> m_values;
477 struct matrix_temporary_type<T,row_major,sparse_bidirectional_iterator_tag> {
482 struct matrix_temporary_type<T,unknown_orientation,sparse_bidirectional_iterator_tag> {
486 template<
class T,
class I>
490 template<
class T,
class I>