compbio
SkylineMatrix.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2009 Guillaume Saupin <guillaume.saupin@cea.fr>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_SKYLINEMATRIX_H
11 #define EIGEN_SKYLINEMATRIX_H
12 
13 #include "SkylineStorage.h"
14 #include "SkylineMatrixBase.h"
15 
16 namespace Eigen {
17 
33 namespace internal {
34 template<typename _Scalar, int _Options>
35 struct traits<SkylineMatrix<_Scalar, _Options> > {
36  typedef _Scalar Scalar;
37  typedef Sparse StorageKind;
38 
39  enum {
40  RowsAtCompileTime = Dynamic,
41  ColsAtCompileTime = Dynamic,
42  MaxRowsAtCompileTime = Dynamic,
43  MaxColsAtCompileTime = Dynamic,
44  Flags = SkylineBit | _Options,
45  CoeffReadCost = NumTraits<Scalar>::ReadCost,
46  };
47 };
48 }
49 
50 template<typename _Scalar, int _Options>
51 class SkylineMatrix
52 : public SkylineMatrixBase<SkylineMatrix<_Scalar, _Options> > {
53 public:
54  EIGEN_SKYLINE_GENERIC_PUBLIC_INTERFACE(SkylineMatrix)
55  EIGEN_SKYLINE_INHERIT_ASSIGNMENT_OPERATOR(SkylineMatrix, +=)
56  EIGEN_SKYLINE_INHERIT_ASSIGNMENT_OPERATOR(SkylineMatrix, -=)
57 
58  using Base::IsRowMajor;
59 
60 protected:
61 
62  typedef SkylineMatrix<Scalar, (Flags&~RowMajorBit) | (IsRowMajor ? RowMajorBit : 0) > TransposedSkylineMatrix;
63 
64  Index m_outerSize;
65  Index m_innerSize;
66 
67 public:
68  Index* m_colStartIndex;
69  Index* m_rowStartIndex;
70  SkylineStorage<Scalar> m_data;
71 
72 public:
73 
74  inline Index rows() const {
75  return IsRowMajor ? m_outerSize : m_innerSize;
76  }
77 
78  inline Index cols() const {
79  return IsRowMajor ? m_innerSize : m_outerSize;
80  }
81 
82  inline Index innerSize() const {
83  return m_innerSize;
84  }
85 
86  inline Index outerSize() const {
87  return m_outerSize;
88  }
89 
90  inline Index upperNonZeros() const {
91  return m_data.upperSize();
92  }
93 
94  inline Index lowerNonZeros() const {
95  return m_data.lowerSize();
96  }
97 
98  inline Index upperNonZeros(Index j) const {
99  return m_colStartIndex[j + 1] - m_colStartIndex[j];
100  }
101 
102  inline Index lowerNonZeros(Index j) const {
103  return m_rowStartIndex[j + 1] - m_rowStartIndex[j];
104  }
105 
106  inline const Scalar* _diagPtr() const {
107  return &m_data.diag(0);
108  }
109 
110  inline Scalar* _diagPtr() {
111  return &m_data.diag(0);
112  }
113 
114  inline const Scalar* _upperPtr() const {
115  return &m_data.upper(0);
116  }
117 
118  inline Scalar* _upperPtr() {
119  return &m_data.upper(0);
120  }
121 
122  inline const Scalar* _lowerPtr() const {
123  return &m_data.lower(0);
124  }
125 
126  inline Scalar* _lowerPtr() {
127  return &m_data.lower(0);
128  }
129 
130  inline const Index* _upperProfilePtr() const {
131  return &m_data.upperProfile(0);
132  }
133 
134  inline Index* _upperProfilePtr() {
135  return &m_data.upperProfile(0);
136  }
137 
138  inline const Index* _lowerProfilePtr() const {
139  return &m_data.lowerProfile(0);
140  }
141 
142  inline Index* _lowerProfilePtr() {
143  return &m_data.lowerProfile(0);
144  }
145 
146  inline Scalar coeff(Index row, Index col) const {
147  const Index outer = IsRowMajor ? row : col;
148  const Index inner = IsRowMajor ? col : row;
149 
150  eigen_assert(outer < outerSize());
151  eigen_assert(inner < innerSize());
152 
153  if (outer == inner)
154  return this->m_data.diag(outer);
155 
156  if (IsRowMajor) {
157  if (inner > outer) //upper matrix
158  {
159  const Index minOuterIndex = inner - m_data.upperProfile(inner);
160  if (outer >= minOuterIndex)
161  return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner)));
162  else
163  return Scalar(0);
164  }
165  if (inner < outer) //lower matrix
166  {
167  const Index minInnerIndex = outer - m_data.lowerProfile(outer);
168  if (inner >= minInnerIndex)
169  return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer)));
170  else
171  return Scalar(0);
172  }
173  return m_data.upper(m_colStartIndex[inner] + outer - inner);
174  } else {
175  if (outer > inner) //upper matrix
176  {
177  const Index maxOuterIndex = inner + m_data.upperProfile(inner);
178  if (outer <= maxOuterIndex)
179  return this->m_data.upper(m_colStartIndex[inner] + (outer - inner));
180  else
181  return Scalar(0);
182  }
183  if (outer < inner) //lower matrix
184  {
185  const Index maxInnerIndex = outer + m_data.lowerProfile(outer);
186 
187  if (inner <= maxInnerIndex)
188  return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer));
189  else
190  return Scalar(0);
191  }
192  }
193  }
194 
195  inline Scalar& coeffRef(Index row, Index col) {
196  const Index outer = IsRowMajor ? row : col;
197  const Index inner = IsRowMajor ? col : row;
198 
199  eigen_assert(outer < outerSize());
200  eigen_assert(inner < innerSize());
201 
202  if (outer == inner)
203  return this->m_data.diag(outer);
204 
205  if (IsRowMajor) {
206  if (col > row) //upper matrix
207  {
208  const Index minOuterIndex = inner - m_data.upperProfile(inner);
209  eigen_assert(outer >= minOuterIndex && "you try to acces a coeff that do not exist in the storage");
210  return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner)));
211  }
212  if (col < row) //lower matrix
213  {
214  const Index minInnerIndex = outer - m_data.lowerProfile(outer);
215  eigen_assert(inner >= minInnerIndex && "you try to acces a coeff that do not exist in the storage");
216  return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer)));
217  }
218  } else {
219  if (outer > inner) //upper matrix
220  {
221  const Index maxOuterIndex = inner + m_data.upperProfile(inner);
222  eigen_assert(outer <= maxOuterIndex && "you try to acces a coeff that do not exist in the storage");
223  return this->m_data.upper(m_colStartIndex[inner] + (outer - inner));
224  }
225  if (outer < inner) //lower matrix
226  {
227  const Index maxInnerIndex = outer + m_data.lowerProfile(outer);
228  eigen_assert(inner <= maxInnerIndex && "you try to acces a coeff that do not exist in the storage");
229  return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer));
230  }
231  }
232  }
233 
234  inline Scalar coeffDiag(Index idx) const {
235  eigen_assert(idx < outerSize());
236  eigen_assert(idx < innerSize());
237  return this->m_data.diag(idx);
238  }
239 
240  inline Scalar coeffLower(Index row, Index col) const {
241  const Index outer = IsRowMajor ? row : col;
242  const Index inner = IsRowMajor ? col : row;
243 
244  eigen_assert(outer < outerSize());
245  eigen_assert(inner < innerSize());
246  eigen_assert(inner != outer);
247 
248  if (IsRowMajor) {
249  const Index minInnerIndex = outer - m_data.lowerProfile(outer);
250  if (inner >= minInnerIndex)
251  return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer)));
252  else
253  return Scalar(0);
254 
255  } else {
256  const Index maxInnerIndex = outer + m_data.lowerProfile(outer);
257  if (inner <= maxInnerIndex)
258  return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer));
259  else
260  return Scalar(0);
261  }
262  }
263 
264  inline Scalar coeffUpper(Index row, Index col) const {
265  const Index outer = IsRowMajor ? row : col;
266  const Index inner = IsRowMajor ? col : row;
267 
268  eigen_assert(outer < outerSize());
269  eigen_assert(inner < innerSize());
270  eigen_assert(inner != outer);
271 
272  if (IsRowMajor) {
273  const Index minOuterIndex = inner - m_data.upperProfile(inner);
274  if (outer >= minOuterIndex)
275  return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner)));
276  else
277  return Scalar(0);
278  } else {
279  const Index maxOuterIndex = inner + m_data.upperProfile(inner);
280  if (outer <= maxOuterIndex)
281  return this->m_data.upper(m_colStartIndex[inner] + (outer - inner));
282  else
283  return Scalar(0);
284  }
285  }
286 
287  inline Scalar& coeffRefDiag(Index idx) {
288  eigen_assert(idx < outerSize());
289  eigen_assert(idx < innerSize());
290  return this->m_data.diag(idx);
291  }
292 
293  inline Scalar& coeffRefLower(Index row, Index col) {
294  const Index outer = IsRowMajor ? row : col;
295  const Index inner = IsRowMajor ? col : row;
296 
297  eigen_assert(outer < outerSize());
298  eigen_assert(inner < innerSize());
299  eigen_assert(inner != outer);
300 
301  if (IsRowMajor) {
302  const Index minInnerIndex = outer - m_data.lowerProfile(outer);
303  eigen_assert(inner >= minInnerIndex && "you try to acces a coeff that do not exist in the storage");
304  return this->m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer)));
305  } else {
306  const Index maxInnerIndex = outer + m_data.lowerProfile(outer);
307  eigen_assert(inner <= maxInnerIndex && "you try to acces a coeff that do not exist in the storage");
308  return this->m_data.lower(m_rowStartIndex[outer] + (inner - outer));
309  }
310  }
311 
312  inline bool coeffExistLower(Index row, Index col) {
313  const Index outer = IsRowMajor ? row : col;
314  const Index inner = IsRowMajor ? col : row;
315 
316  eigen_assert(outer < outerSize());
317  eigen_assert(inner < innerSize());
318  eigen_assert(inner != outer);
319 
320  if (IsRowMajor) {
321  const Index minInnerIndex = outer - m_data.lowerProfile(outer);
322  return inner >= minInnerIndex;
323  } else {
324  const Index maxInnerIndex = outer + m_data.lowerProfile(outer);
325  return inner <= maxInnerIndex;
326  }
327  }
328 
329  inline Scalar& coeffRefUpper(Index row, Index col) {
330  const Index outer = IsRowMajor ? row : col;
331  const Index inner = IsRowMajor ? col : row;
332 
333  eigen_assert(outer < outerSize());
334  eigen_assert(inner < innerSize());
335  eigen_assert(inner != outer);
336 
337  if (IsRowMajor) {
338  const Index minOuterIndex = inner - m_data.upperProfile(inner);
339  eigen_assert(outer >= minOuterIndex && "you try to acces a coeff that do not exist in the storage");
340  return this->m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner)));
341  } else {
342  const Index maxOuterIndex = inner + m_data.upperProfile(inner);
343  eigen_assert(outer <= maxOuterIndex && "you try to acces a coeff that do not exist in the storage");
344  return this->m_data.upper(m_colStartIndex[inner] + (outer - inner));
345  }
346  }
347 
348  inline bool coeffExistUpper(Index row, Index col) {
349  const Index outer = IsRowMajor ? row : col;
350  const Index inner = IsRowMajor ? col : row;
351 
352  eigen_assert(outer < outerSize());
353  eigen_assert(inner < innerSize());
354  eigen_assert(inner != outer);
355 
356  if (IsRowMajor) {
357  const Index minOuterIndex = inner - m_data.upperProfile(inner);
358  return outer >= minOuterIndex;
359  } else {
360  const Index maxOuterIndex = inner + m_data.upperProfile(inner);
361  return outer <= maxOuterIndex;
362  }
363  }
364 
365 
366 protected:
367 
368 public:
369  class InnerUpperIterator;
370  class InnerLowerIterator;
371 
372  class OuterUpperIterator;
373  class OuterLowerIterator;
374 
376  inline void setZero() {
377  m_data.clear();
378  memset(m_colStartIndex, 0, (m_outerSize + 1) * sizeof (Index));
379  memset(m_rowStartIndex, 0, (m_outerSize + 1) * sizeof (Index));
380  }
381 
383  inline Index nonZeros() const {
384  return m_data.diagSize() + m_data.upperSize() + m_data.lowerSize();
385  }
386 
388  inline void reserve(Index reserveSize, Index reserveUpperSize, Index reserveLowerSize) {
389  m_data.reserve(reserveSize, reserveUpperSize, reserveLowerSize);
390  }
391 
400  EIGEN_DONT_INLINE Scalar & insert(Index row, Index col) {
401  const Index outer = IsRowMajor ? row : col;
402  const Index inner = IsRowMajor ? col : row;
403 
404  eigen_assert(outer < outerSize());
405  eigen_assert(inner < innerSize());
406 
407  if (outer == inner)
408  return m_data.diag(col);
409 
410  if (IsRowMajor) {
411  if (outer < inner) //upper matrix
412  {
413  Index minOuterIndex = 0;
414  minOuterIndex = inner - m_data.upperProfile(inner);
415 
416  if (outer < minOuterIndex) //The value does not yet exist
417  {
418  const Index previousProfile = m_data.upperProfile(inner);
419 
420  m_data.upperProfile(inner) = inner - outer;
421 
422 
423  const Index bandIncrement = m_data.upperProfile(inner) - previousProfile;
424  //shift data stored after this new one
425  const Index stop = m_colStartIndex[cols()];
426  const Index start = m_colStartIndex[inner];
427 
428 
429  for (Index innerIdx = stop; innerIdx >= start; innerIdx--) {
430  m_data.upper(innerIdx + bandIncrement) = m_data.upper(innerIdx);
431  }
432 
433  for (Index innerIdx = cols(); innerIdx > inner; innerIdx--) {
434  m_colStartIndex[innerIdx] += bandIncrement;
435  }
436 
437  //zeros new data
438  memset(this->_upperPtr() + start, 0, (bandIncrement - 1) * sizeof (Scalar));
439 
440  return m_data.upper(m_colStartIndex[inner]);
441  } else {
442  return m_data.upper(m_colStartIndex[inner] + outer - (inner - m_data.upperProfile(inner)));
443  }
444  }
445 
446  if (outer > inner) //lower matrix
447  {
448  const Index minInnerIndex = outer - m_data.lowerProfile(outer);
449  if (inner < minInnerIndex) //The value does not yet exist
450  {
451  const Index previousProfile = m_data.lowerProfile(outer);
452  m_data.lowerProfile(outer) = outer - inner;
453 
454  const Index bandIncrement = m_data.lowerProfile(outer) - previousProfile;
455  //shift data stored after this new one
456  const Index stop = m_rowStartIndex[rows()];
457  const Index start = m_rowStartIndex[outer];
458 
459 
460  for (Index innerIdx = stop; innerIdx >= start; innerIdx--) {
461  m_data.lower(innerIdx + bandIncrement) = m_data.lower(innerIdx);
462  }
463 
464  for (Index innerIdx = rows(); innerIdx > outer; innerIdx--) {
465  m_rowStartIndex[innerIdx] += bandIncrement;
466  }
467 
468  //zeros new data
469  memset(this->_lowerPtr() + start, 0, (bandIncrement - 1) * sizeof (Scalar));
470  return m_data.lower(m_rowStartIndex[outer]);
471  } else {
472  return m_data.lower(m_rowStartIndex[outer] + inner - (outer - m_data.lowerProfile(outer)));
473  }
474  }
475  } else {
476  if (outer > inner) //upper matrix
477  {
478  const Index maxOuterIndex = inner + m_data.upperProfile(inner);
479  if (outer > maxOuterIndex) //The value does not yet exist
480  {
481  const Index previousProfile = m_data.upperProfile(inner);
482  m_data.upperProfile(inner) = outer - inner;
483 
484  const Index bandIncrement = m_data.upperProfile(inner) - previousProfile;
485  //shift data stored after this new one
486  const Index stop = m_rowStartIndex[rows()];
487  const Index start = m_rowStartIndex[inner + 1];
488 
489  for (Index innerIdx = stop; innerIdx >= start; innerIdx--) {
490  m_data.upper(innerIdx + bandIncrement) = m_data.upper(innerIdx);
491  }
492 
493  for (Index innerIdx = inner + 1; innerIdx < outerSize() + 1; innerIdx++) {
494  m_rowStartIndex[innerIdx] += bandIncrement;
495  }
496  memset(this->_upperPtr() + m_rowStartIndex[inner] + previousProfile + 1, 0, (bandIncrement - 1) * sizeof (Scalar));
497  return m_data.upper(m_rowStartIndex[inner] + m_data.upperProfile(inner));
498  } else {
499  return m_data.upper(m_rowStartIndex[inner] + (outer - inner));
500  }
501  }
502 
503  if (outer < inner) //lower matrix
504  {
505  const Index maxInnerIndex = outer + m_data.lowerProfile(outer);
506  if (inner > maxInnerIndex) //The value does not yet exist
507  {
508  const Index previousProfile = m_data.lowerProfile(outer);
509  m_data.lowerProfile(outer) = inner - outer;
510 
511  const Index bandIncrement = m_data.lowerProfile(outer) - previousProfile;
512  //shift data stored after this new one
513  const Index stop = m_colStartIndex[cols()];
514  const Index start = m_colStartIndex[outer + 1];
515 
516  for (Index innerIdx = stop; innerIdx >= start; innerIdx--) {
517  m_data.lower(innerIdx + bandIncrement) = m_data.lower(innerIdx);
518  }
519 
520  for (Index innerIdx = outer + 1; innerIdx < outerSize() + 1; innerIdx++) {
521  m_colStartIndex[innerIdx] += bandIncrement;
522  }
523  memset(this->_lowerPtr() + m_colStartIndex[outer] + previousProfile + 1, 0, (bandIncrement - 1) * sizeof (Scalar));
524  return m_data.lower(m_colStartIndex[outer] + m_data.lowerProfile(outer));
525  } else {
526  return m_data.lower(m_colStartIndex[outer] + (inner - outer));
527  }
528  }
529  }
530  }
531 
534  inline void finalize() {
535  if (IsRowMajor) {
536  if (rows() > cols())
537  m_data.resize(cols(), cols(), rows(), m_colStartIndex[cols()] + 1, m_rowStartIndex[rows()] + 1);
538  else
539  m_data.resize(rows(), cols(), rows(), m_colStartIndex[cols()] + 1, m_rowStartIndex[rows()] + 1);
540 
541  // eigen_assert(rows() == cols() && "memory reorganisatrion only works with suare matrix");
542  //
543  // Scalar* newArray = new Scalar[m_colStartIndex[cols()] + 1 + m_rowStartIndex[rows()] + 1];
544  // Index dataIdx = 0;
545  // for (Index row = 0; row < rows(); row++) {
546  //
547  // const Index nbLowerElts = m_rowStartIndex[row + 1] - m_rowStartIndex[row];
548  // // std::cout << "nbLowerElts" << nbLowerElts << std::endl;
549  // memcpy(newArray + dataIdx, m_data.m_lower + m_rowStartIndex[row], nbLowerElts * sizeof (Scalar));
550  // m_rowStartIndex[row] = dataIdx;
551  // dataIdx += nbLowerElts;
552  //
553  // const Index nbUpperElts = m_colStartIndex[row + 1] - m_colStartIndex[row];
554  // memcpy(newArray + dataIdx, m_data.m_upper + m_colStartIndex[row], nbUpperElts * sizeof (Scalar));
555  // m_colStartIndex[row] = dataIdx;
556  // dataIdx += nbUpperElts;
557  //
558  //
559  // }
560  // //todo : don't access m_data profile directly : add an accessor from SkylineMatrix
561  // m_rowStartIndex[rows()] = m_rowStartIndex[rows()-1] + m_data.lowerProfile(rows()-1);
562  // m_colStartIndex[cols()] = m_colStartIndex[cols()-1] + m_data.upperProfile(cols()-1);
563  //
564  // delete[] m_data.m_lower;
565  // delete[] m_data.m_upper;
566  //
567  // m_data.m_lower = newArray;
568  // m_data.m_upper = newArray;
569  } else {
570  if (rows() > cols())
571  m_data.resize(cols(), rows(), cols(), m_rowStartIndex[cols()] + 1, m_colStartIndex[cols()] + 1);
572  else
573  m_data.resize(rows(), rows(), cols(), m_rowStartIndex[rows()] + 1, m_colStartIndex[rows()] + 1);
574  }
575  }
576 
577  inline void squeeze() {
578  finalize();
579  m_data.squeeze();
580  }
581 
582  void prune(Scalar reference, RealScalar epsilon = dummy_precision<RealScalar > ()) {
583  //TODO
584  }
585 
589  void resize(size_t rows, size_t cols) {
590  const Index diagSize = rows > cols ? cols : rows;
591  m_innerSize = IsRowMajor ? cols : rows;
592 
593  eigen_assert(rows == cols && "Skyline matrix must be square matrix");
594 
595  if (diagSize % 2) { // diagSize is odd
596  const Index k = (diagSize - 1) / 2;
597 
598  m_data.resize(diagSize, IsRowMajor ? cols : rows, IsRowMajor ? rows : cols,
599  2 * k * k + k + 1,
600  2 * k * k + k + 1);
601 
602  } else // diagSize is even
603  {
604  const Index k = diagSize / 2;
605  m_data.resize(diagSize, IsRowMajor ? cols : rows, IsRowMajor ? rows : cols,
606  2 * k * k - k + 1,
607  2 * k * k - k + 1);
608  }
609 
610  if (m_colStartIndex && m_rowStartIndex) {
611  delete[] m_colStartIndex;
612  delete[] m_rowStartIndex;
613  }
614  m_colStartIndex = new Index [cols + 1];
615  m_rowStartIndex = new Index [rows + 1];
616  m_outerSize = diagSize;
617 
618  m_data.reset();
619  m_data.clear();
620 
621  m_outerSize = diagSize;
622  memset(m_colStartIndex, 0, (cols + 1) * sizeof (Index));
623  memset(m_rowStartIndex, 0, (rows + 1) * sizeof (Index));
624  }
625 
626  void resizeNonZeros(Index size) {
627  m_data.resize(size);
628  }
629 
630  inline SkylineMatrix()
631  : m_outerSize(-1), m_innerSize(0), m_colStartIndex(0), m_rowStartIndex(0) {
632  resize(0, 0);
633  }
634 
635  inline SkylineMatrix(size_t rows, size_t cols)
636  : m_outerSize(0), m_innerSize(0), m_colStartIndex(0), m_rowStartIndex(0) {
637  resize(rows, cols);
638  }
639 
640  template<typename OtherDerived>
641  inline SkylineMatrix(const SkylineMatrixBase<OtherDerived>& other)
642  : m_outerSize(0), m_innerSize(0), m_colStartIndex(0), m_rowStartIndex(0) {
643  *this = other.derived();
644  }
645 
646  inline SkylineMatrix(const SkylineMatrix & other)
647  : Base(), m_outerSize(0), m_innerSize(0), m_colStartIndex(0), m_rowStartIndex(0) {
648  *this = other.derived();
649  }
650 
651  inline void swap(SkylineMatrix & other) {
652  //EIGEN_DBG_SKYLINE(std::cout << "SkylineMatrix:: swap\n");
653  std::swap(m_colStartIndex, other.m_colStartIndex);
654  std::swap(m_rowStartIndex, other.m_rowStartIndex);
655  std::swap(m_innerSize, other.m_innerSize);
656  std::swap(m_outerSize, other.m_outerSize);
657  m_data.swap(other.m_data);
658  }
659 
660  inline SkylineMatrix & operator=(const SkylineMatrix & other) {
661  std::cout << "SkylineMatrix& operator=(const SkylineMatrix& other)\n";
662  if (other.isRValue()) {
663  swap(other.const_cast_derived());
664  } else {
665  resize(other.rows(), other.cols());
666  memcpy(m_colStartIndex, other.m_colStartIndex, (m_outerSize + 1) * sizeof (Index));
667  memcpy(m_rowStartIndex, other.m_rowStartIndex, (m_outerSize + 1) * sizeof (Index));
668  m_data = other.m_data;
669  }
670  return *this;
671  }
672 
673  template<typename OtherDerived>
674  inline SkylineMatrix & operator=(const SkylineMatrixBase<OtherDerived>& other) {
675  const bool needToTranspose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit);
676  if (needToTranspose) {
677  // TODO
678  // return *this;
679  } else {
680  // there is no special optimization
681  return SkylineMatrixBase<SkylineMatrix>::operator=(other.derived());
682  }
683  }
684 
685  friend std::ostream & operator <<(std::ostream & s, const SkylineMatrix & m) {
686 
687  EIGEN_DBG_SKYLINE(
688  std::cout << "upper elements : " << std::endl;
689  for (Index i = 0; i < m.m_data.upperSize(); i++)
690  std::cout << m.m_data.upper(i) << "\t";
691  std::cout << std::endl;
692  std::cout << "upper profile : " << std::endl;
693  for (Index i = 0; i < m.m_data.upperProfileSize(); i++)
694  std::cout << m.m_data.upperProfile(i) << "\t";
695  std::cout << std::endl;
696  std::cout << "lower startIdx : " << std::endl;
697  for (Index i = 0; i < m.m_data.upperProfileSize(); i++)
698  std::cout << (IsRowMajor ? m.m_colStartIndex[i] : m.m_rowStartIndex[i]) << "\t";
699  std::cout << std::endl;
700 
701 
702  std::cout << "lower elements : " << std::endl;
703  for (Index i = 0; i < m.m_data.lowerSize(); i++)
704  std::cout << m.m_data.lower(i) << "\t";
705  std::cout << std::endl;
706  std::cout << "lower profile : " << std::endl;
707  for (Index i = 0; i < m.m_data.lowerProfileSize(); i++)
708  std::cout << m.m_data.lowerProfile(i) << "\t";
709  std::cout << std::endl;
710  std::cout << "lower startIdx : " << std::endl;
711  for (Index i = 0; i < m.m_data.lowerProfileSize(); i++)
712  std::cout << (IsRowMajor ? m.m_rowStartIndex[i] : m.m_colStartIndex[i]) << "\t";
713  std::cout << std::endl;
714  );
715  for (Index rowIdx = 0; rowIdx < m.rows(); rowIdx++) {
716  for (Index colIdx = 0; colIdx < m.cols(); colIdx++) {
717  s << m.coeff(rowIdx, colIdx) << "\t";
718  }
719  s << std::endl;
720  }
721  return s;
722  }
723 
725  inline ~SkylineMatrix() {
726  delete[] m_colStartIndex;
727  delete[] m_rowStartIndex;
728  }
729 
731  Scalar sum() const;
732 };
733 
734 template<typename Scalar, int _Options>
735 class SkylineMatrix<Scalar, _Options>::InnerUpperIterator {
736 public:
737 
738  InnerUpperIterator(const SkylineMatrix& mat, Index outer)
739  : m_matrix(mat), m_outer(outer),
740  m_id(_Options == RowMajor ? mat.m_colStartIndex[outer] : mat.m_rowStartIndex[outer] + 1),
741  m_start(m_id),
742  m_end(_Options == RowMajor ? mat.m_colStartIndex[outer + 1] : mat.m_rowStartIndex[outer + 1] + 1) {
743  }
744 
745  inline InnerUpperIterator & operator++() {
746  m_id++;
747  return *this;
748  }
749 
750  inline InnerUpperIterator & operator+=(Index shift) {
751  m_id += shift;
752  return *this;
753  }
754 
755  inline Scalar value() const {
756  return m_matrix.m_data.upper(m_id);
757  }
758 
759  inline Scalar* valuePtr() {
760  return const_cast<Scalar*> (&(m_matrix.m_data.upper(m_id)));
761  }
762 
763  inline Scalar& valueRef() {
764  return const_cast<Scalar&> (m_matrix.m_data.upper(m_id));
765  }
766 
767  inline Index index() const {
768  return IsRowMajor ? m_outer - m_matrix.m_data.upperProfile(m_outer) + (m_id - m_start) :
769  m_outer + (m_id - m_start) + 1;
770  }
771 
772  inline Index row() const {
773  return IsRowMajor ? index() : m_outer;
774  }
775 
776  inline Index col() const {
777  return IsRowMajor ? m_outer : index();
778  }
779 
780  inline size_t size() const {
781  return m_matrix.m_data.upperProfile(m_outer);
782  }
783 
784  inline operator bool() const {
785  return (m_id < m_end) && (m_id >= m_start);
786  }
787 
788 protected:
789  const SkylineMatrix& m_matrix;
790  const Index m_outer;
791  Index m_id;
792  const Index m_start;
793  const Index m_end;
794 };
795 
796 template<typename Scalar, int _Options>
797 class SkylineMatrix<Scalar, _Options>::InnerLowerIterator {
798 public:
799 
800  InnerLowerIterator(const SkylineMatrix& mat, Index outer)
801  : m_matrix(mat),
802  m_outer(outer),
803  m_id(_Options == RowMajor ? mat.m_rowStartIndex[outer] : mat.m_colStartIndex[outer] + 1),
804  m_start(m_id),
805  m_end(_Options == RowMajor ? mat.m_rowStartIndex[outer + 1] : mat.m_colStartIndex[outer + 1] + 1) {
806  }
807 
808  inline InnerLowerIterator & operator++() {
809  m_id++;
810  return *this;
811  }
812 
813  inline InnerLowerIterator & operator+=(Index shift) {
814  m_id += shift;
815  return *this;
816  }
817 
818  inline Scalar value() const {
819  return m_matrix.m_data.lower(m_id);
820  }
821 
822  inline Scalar* valuePtr() {
823  return const_cast<Scalar*> (&(m_matrix.m_data.lower(m_id)));
824  }
825 
826  inline Scalar& valueRef() {
827  return const_cast<Scalar&> (m_matrix.m_data.lower(m_id));
828  }
829 
830  inline Index index() const {
831  return IsRowMajor ? m_outer - m_matrix.m_data.lowerProfile(m_outer) + (m_id - m_start) :
832  m_outer + (m_id - m_start) + 1;
833  ;
834  }
835 
836  inline Index row() const {
837  return IsRowMajor ? m_outer : index();
838  }
839 
840  inline Index col() const {
841  return IsRowMajor ? index() : m_outer;
842  }
843 
844  inline size_t size() const {
845  return m_matrix.m_data.lowerProfile(m_outer);
846  }
847 
848  inline operator bool() const {
849  return (m_id < m_end) && (m_id >= m_start);
850  }
851 
852 protected:
853  const SkylineMatrix& m_matrix;
854  const Index m_outer;
855  Index m_id;
856  const Index m_start;
857  const Index m_end;
858 };
859 
860 } // end namespace Eigen
861 
862 #endif // EIGEN_SkylineMatrix_H
Eigen::SkylineStorage
Stores a skyline set of values in three structures : The diagonal elements The upper elements The low...
Definition: SkylineStorage.h:22
Eigen::SkylineMatrix::setZero
void setZero()
Removes all non zeros.
Definition: SkylineMatrix.h:376
Eigen
Namespace containing all symbols from the Eigen library.
Definition: bench_norm.cpp:85
Eigen::SkylineMatrix::nonZeros
Index nonZeros() const
Definition: SkylineMatrix.h:383
Eigen::NumTraits
Holds information about the various numeric (i.e.
Definition: NumTraits.h:150
Eigen::SkylineMatrix::resize
void resize(size_t rows, size_t cols)
Resizes the matrix to a rows x cols matrix and initializes it to zero.
Definition: SkylineMatrix.h:589
Eigen::EigenBase::Index
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:37
Eigen::RowMajorBit
const unsigned int RowMajorBit
for a matrix, this means that the storage order is row-major.
Definition: Constants.h:61
Eigen::SkylineMatrix::reserve
void reserve(Index reserveSize, Index reserveUpperSize, Index reserveLowerSize)
Preallocates reserveSize non zeros.
Definition: SkylineMatrix.h:388
Eigen::SkylineMatrix::insert
EIGEN_DONT_INLINE Scalar & insert(Index row, Index col)
Definition: SkylineMatrix.h:400
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
Eigen::Sparse
The type used to identify a general sparse storage.
Definition: Constants.h:494
Eigen::SkylineMatrix::finalize
void finalize()
Must be called after inserting a set of non zero entries.
Definition: SkylineMatrix.h:534
internal
Definition: BandTriangularSolver.h:13
Eigen::RowMajor
Storage order is row major (see TopicStorageOrders).
Definition: Constants.h:322
Eigen::SkylineMatrixBase
Base class of any skyline matrices or skyline expressions.
Definition: SkylineMatrixBase.h:26
Eigen::SkylineMatrix
The main skyline matrix class.
Definition: SkylineMatrix.h:51
Eigen::SkylineMatrix::InnerUpperIterator
Definition: SkylineMatrix.h:735
Eigen::Dynamic
const int Dynamic
This value means that a positive quantity (e.g., a size) is not known at compile-time, and that instead the value is stored in some runtime variable.
Definition: Constants.h:21
Eigen::SkylineMatrix::~SkylineMatrix
~SkylineMatrix()
Destructor.
Definition: SkylineMatrix.h:725
Eigen::SkylineMatrixBase< SkylineMatrix< _Scalar, _Options > >::RealScalar
NumTraits< Scalar >::Real RealScalar
This is the "real scalar" type; if the Scalar type is already real numbers (e.g.
Definition: SkylineMatrixBase.h:85
Eigen::internal::traits
Definition: ForwardDeclarations.h:17
Eigen::SkylineMatrix::InnerLowerIterator
Definition: SkylineMatrix.h:797
Generated by   doxygen 1.8.12