BlockSparseMatrix.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2013 Desire Nuentsa <desire.nuentsa_wakam@inria.fr>
5 // Copyright (C) 2013 Gael Guennebaud <gael.guennebaud@inria.fr>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11 #ifndef EIGEN_SPARSEBLOCKMATRIX_H
12 #define EIGEN_SPARSEBLOCKMATRIX_H
13 
14 namespace Eigen {
54 template<typename _Scalar, int _BlockAtCompileTime=Dynamic, int _Options=ColMajor, typename _StorageIndex=int> class BlockSparseMatrix;
55 
56 template<typename BlockSparseMatrixT> class BlockSparseMatrixView;
57 
58 namespace internal {
59 template<typename _Scalar, int _BlockAtCompileTime, int _Options, typename _Index>
60 struct traits<BlockSparseMatrix<_Scalar,_BlockAtCompileTime,_Options, _Index> >
61 {
62  typedef _Scalar Scalar;
63  typedef _Index Index;
64  typedef Sparse StorageKind; // FIXME Where is it used ??
65  typedef MatrixXpr XprKind;
66  enum {
67  RowsAtCompileTime = Dynamic,
68  ColsAtCompileTime = Dynamic,
69  MaxRowsAtCompileTime = Dynamic,
70  MaxColsAtCompileTime = Dynamic,
71  BlockSize = _BlockAtCompileTime,
72  Flags = _Options | NestByRefBit | LvalueBit,
73  CoeffReadCost = NumTraits<Scalar>::ReadCost,
74  SupportedAccessPatterns = InnerRandomAccessPattern
75  };
76 };
77 template<typename BlockSparseMatrixT>
78 struct traits<BlockSparseMatrixView<BlockSparseMatrixT> >
79 {
82 
83 };
84 
85 // Function object to sort a triplet list
86 template<typename Iterator, bool IsColMajor>
88 {
89  typedef typename Iterator::value_type Triplet;
90  bool operator()(const Triplet& a, const Triplet& b)
91  { if(IsColMajor)
92  return ((a.col() == b.col() && a.row() < b.row()) || (a.col() < b.col()));
93  else
94  return ((a.row() == b.row() && a.col() < b.col()) || (a.row() < b.row()));
95  }
96 };
97 } // end namespace internal
98 
99 
100 /* Proxy to view the block sparse matrix as a regular sparse matrix */
101 template<typename BlockSparseMatrixT>
102 class BlockSparseMatrixView : public SparseMatrixBase<BlockSparseMatrixT>
103 {
104  public:
108  typedef BlockSparseMatrixT Nested;
109  enum {
110  Flags = BlockSparseMatrixT::Options,
111  Options = BlockSparseMatrixT::Options,
112  RowsAtCompileTime = BlockSparseMatrixT::RowsAtCompileTime,
113  ColsAtCompileTime = BlockSparseMatrixT::ColsAtCompileTime,
114  MaxColsAtCompileTime = BlockSparseMatrixT::MaxColsAtCompileTime,
115  MaxRowsAtCompileTime = BlockSparseMatrixT::MaxRowsAtCompileTime
116  };
117  public:
118  BlockSparseMatrixView(const BlockSparseMatrixT& spblockmat)
119  : m_spblockmat(spblockmat)
120  {}
121 
122  Index outerSize() const
123  {
124  return (Flags&RowMajorBit) == 1 ? this->rows() : this->cols();
125  }
126  Index cols() const
127  {
128  return m_spblockmat.blockCols();
129  }
130  Index rows() const
131  {
132  return m_spblockmat.blockRows();
133  }
134  Scalar coeff(Index row, Index col)
135  {
136  return m_spblockmat.coeff(row, col);
137  }
138  Scalar coeffRef(Index row, Index col)
139  {
140  return m_spblockmat.coeffRef(row, col);
141  }
142  // Wrapper to iterate over all blocks
143  class InnerIterator : public BlockSparseMatrixT::BlockInnerIterator
144  {
145  public:
147  : BlockSparseMatrixT::BlockInnerIterator(mat.m_spblockmat, outer)
148  {}
149 
150  };
151 
152  protected:
153  const BlockSparseMatrixT& m_spblockmat;
154 };
155 
156 // Proxy to view a regular vector as a block vector
157 template<typename BlockSparseMatrixT, typename VectorType>
159 {
160  public:
161  enum {
162  BlockSize = BlockSparseMatrixT::BlockSize,
163  ColsAtCompileTime = VectorType::ColsAtCompileTime,
164  RowsAtCompileTime = VectorType::RowsAtCompileTime,
165  Flags = VectorType::Flags
166  };
169  public:
170  BlockVectorView(const BlockSparseMatrixT& spblockmat, const VectorType& vec)
171  : m_spblockmat(spblockmat),m_vec(vec)
172  { }
173  inline Index cols() const
174  {
175  return m_vec.cols();
176  }
177  inline Index size() const
178  {
179  return m_spblockmat.blockRows();
180  }
181  inline Scalar coeff(Index bi) const
182  {
183  Index startRow = m_spblockmat.blockRowsIndex(bi);
184  Index rowSize = m_spblockmat.blockRowsIndex(bi+1) - startRow;
185  return m_vec.middleRows(startRow, rowSize);
186  }
187  inline Scalar coeff(Index bi, Index j) const
188  {
189  Index startRow = m_spblockmat.blockRowsIndex(bi);
190  Index rowSize = m_spblockmat.blockRowsIndex(bi+1) - startRow;
191  return m_vec.block(startRow, j, rowSize, 1);
192  }
193  protected:
194  const BlockSparseMatrixT& m_spblockmat;
196 };
197 
198 template<typename VectorType, typename Index> class BlockVectorReturn;
199 
200 
201 // Proxy to view a regular vector as a block vector
202 template<typename BlockSparseMatrixT, typename VectorType>
203 class BlockVectorReturn
204 {
205  public:
206  enum {
207  ColsAtCompileTime = VectorType::ColsAtCompileTime,
208  RowsAtCompileTime = VectorType::RowsAtCompileTime,
209  Flags = VectorType::Flags
210  };
213  public:
214  BlockVectorReturn(const BlockSparseMatrixT& spblockmat, VectorType& vec)
215  : m_spblockmat(spblockmat),m_vec(vec)
216  { }
217  inline Index size() const
218  {
219  return m_spblockmat.blockRows();
220  }
221  inline Scalar coeffRef(Index bi)
222  {
223  Index startRow = m_spblockmat.blockRowsIndex(bi);
224  Index rowSize = m_spblockmat.blockRowsIndex(bi+1) - startRow;
225  return m_vec.middleRows(startRow, rowSize);
226  }
227  inline Scalar coeffRef(Index bi, Index j)
228  {
229  Index startRow = m_spblockmat.blockRowsIndex(bi);
230  Index rowSize = m_spblockmat.blockRowsIndex(bi+1) - startRow;
231  return m_vec.block(startRow, j, rowSize, 1);
232  }
233 
234  protected:
235  const BlockSparseMatrixT& m_spblockmat;
237 };
238 
239 // Block version of the sparse dense product
240 template<typename Lhs, typename Rhs>
242 
243 namespace internal {
244 
245 template<typename BlockSparseMatrixT, typename VecType>
246 struct traits<BlockSparseTimeDenseProduct<BlockSparseMatrixT, VecType> >
247 {
252  enum {
253  RowsAtCompileTime = Dynamic,
254  ColsAtCompileTime = Dynamic,
255  MaxRowsAtCompileTime = Dynamic,
256  MaxColsAtCompileTime = Dynamic,
257  Flags = 0,
259  };
260 };
261 } // end namespace internal
262 
263 template<typename Lhs, typename Rhs>
265  : public ProductBase<BlockSparseTimeDenseProduct<Lhs,Rhs>, Lhs, Rhs>
266 {
267  public:
268  EIGEN_PRODUCT_PUBLIC_INTERFACE(BlockSparseTimeDenseProduct)
269 
270  BlockSparseTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
271  {}
272 
273  template<typename Dest> void scaleAndAddTo(Dest& dest, const typename Rhs::Scalar& alpha) const
274  {
275  BlockVectorReturn<Lhs,Dest> tmpDest(m_lhs, dest);
277  }
278 
279  private:
281 };
282 
283 template<typename _Scalar, int _BlockAtCompileTime, int _Options, typename _StorageIndex>
284 class BlockSparseMatrix : public SparseMatrixBase<BlockSparseMatrix<_Scalar,_BlockAtCompileTime, _Options,_StorageIndex> >
285 {
286  public:
287  typedef _Scalar Scalar;
289  typedef _StorageIndex StorageIndex;
291 
292  enum {
293  Options = _Options,
294  Flags = Options,
295  BlockSize=_BlockAtCompileTime,
296  RowsAtCompileTime = Dynamic,
297  ColsAtCompileTime = Dynamic,
298  MaxRowsAtCompileTime = Dynamic,
299  MaxColsAtCompileTime = Dynamic,
300  IsVectorAtCompileTime = 0,
301  IsColMajor = Flags&RowMajorBit ? 0 : 1
302  };
307  public:
308  // Default constructor
310  : m_innerBSize(0),m_outerBSize(0),m_innerOffset(0),m_outerOffset(0),
311  m_nonzerosblocks(0),m_values(0),m_blockPtr(0),m_indices(0),
312  m_outerIndex(0),m_blockSize(BlockSize)
313  { }
314 
315 
321  : m_innerBSize(IsColMajor ? brow : bcol),
322  m_outerBSize(IsColMajor ? bcol : brow),
323  m_innerOffset(0),m_outerOffset(0),m_nonzerosblocks(0),
324  m_values(0),m_blockPtr(0),m_indices(0),
325  m_outerIndex(0),m_blockSize(BlockSize)
326  { }
327 
332  : m_innerBSize(other.m_innerBSize),m_outerBSize(other.m_outerBSize),
333  m_nonzerosblocks(other.m_nonzerosblocks),m_nonzeros(other.m_nonzeros),
334  m_blockPtr(0),m_blockSize(other.m_blockSize)
335  {
336  // should we allow copying between variable-size blocks and fixed-size blocks ??
337  eigen_assert(m_blockSize == BlockSize && " CAN NOT COPY BETWEEN FIXED-SIZE AND VARIABLE-SIZE BLOCKS");
338 
339  std::copy(other.m_innerOffset, other.m_innerOffset+m_innerBSize+1, m_innerOffset);
340  std::copy(other.m_outerOffset, other.m_outerOffset+m_outerBSize+1, m_outerOffset);
341  std::copy(other.m_values, other.m_values+m_nonzeros, m_values);
342 
343  if(m_blockSize != Dynamic)
344  std::copy(other.m_blockPtr, other.m_blockPtr+m_nonzerosblocks, m_blockPtr);
345 
346  std::copy(other.m_indices, other.m_indices+m_nonzerosblocks, m_indices);
347  std::copy(other.m_outerIndex, other.m_outerIndex+m_outerBSize, m_outerIndex);
348  }
349 
351  {
352  std::swap(first.m_innerBSize, second.m_innerBSize);
353  std::swap(first.m_outerBSize, second.m_outerBSize);
354  std::swap(first.m_innerOffset, second.m_innerOffset);
355  std::swap(first.m_outerOffset, second.m_outerOffset);
357  std::swap(first.m_nonzeros, second.m_nonzeros);
358  std::swap(first.m_values, second.m_values);
359  std::swap(first.m_blockPtr, second.m_blockPtr);
360  std::swap(first.m_indices, second.m_indices);
361  std::swap(first.m_outerIndex, second.m_outerIndex);
362  std::swap(first.m_BlockSize, second.m_blockSize);
363  }
364 
366  {
367  //Copy-and-swap paradigm ... avoid leaked data if thrown
368  swap(*this, other);
369  return *this;
370  }
371 
372  // Destructor
374  {
375  delete[] m_outerIndex;
376  delete[] m_innerOffset;
377  delete[] m_outerOffset;
378  delete[] m_indices;
379  delete[] m_blockPtr;
380  delete[] m_values;
381  }
382 
383 
388  template<typename MatrixType>
389  inline BlockSparseMatrix(const MatrixType& spmat) : m_blockSize(BlockSize)
390  {
391  EIGEN_STATIC_ASSERT((m_blockSize != Dynamic), THIS_METHOD_IS_ONLY_FOR_FIXED_SIZE);
392 
393  *this = spmat;
394  }
395 
404  template<typename MatrixType>
405  inline BlockSparseMatrix& operator=(const MatrixType& spmat)
406  {
407  eigen_assert((m_innerBSize != 0 && m_outerBSize != 0)
408  && "Trying to assign to a zero-size matrix, call resize() first");
409  eigen_assert(((MatrixType::Options&RowMajorBit) != IsColMajor) && "Wrong storage order");
411  MatrixPatternType blockPattern(blockRows(), blockCols());
412  m_nonzeros = 0;
413 
414  // First, compute the number of nonzero blocks and their locations
415  for(StorageIndex bj = 0; bj < m_outerBSize; ++bj)
416  {
417  // Browse each outer block and compute the structure
418  std::vector<bool> nzblocksFlag(m_innerBSize,false); // Record the existing blocks
419  blockPattern.startVec(bj);
420  for(StorageIndex j = blockOuterIndex(bj); j < blockOuterIndex(bj+1); ++j)
421  {
422  typename MatrixType::InnerIterator it_spmat(spmat, j);
423  for(; it_spmat; ++it_spmat)
424  {
425  StorageIndex bi = innerToBlock(it_spmat.index()); // Index of the current nonzero block
426  if(!nzblocksFlag[bi])
427  {
428  // Save the index of this nonzero block
429  nzblocksFlag[bi] = true;
430  blockPattern.insertBackByOuterInnerUnordered(bj, bi) = true;
431  // Compute the total number of nonzeros (including explicit zeros in blocks)
432  m_nonzeros += blockOuterSize(bj) * blockInnerSize(bi);
433  }
434  }
435  } // end current outer block
436  }
437  blockPattern.finalize();
438 
439  // Allocate the internal arrays
440  setBlockStructure(blockPattern);
441 
442  for(StorageIndex nz = 0; nz < m_nonzeros; ++nz) m_values[nz] = Scalar(0);
443  for(StorageIndex bj = 0; bj < m_outerBSize; ++bj)
444  {
445  // Now copy the values
446  for(StorageIndex j = blockOuterIndex(bj); j < blockOuterIndex(bj+1); ++j)
447  {
448  // Browse the outer block column by column (for column-major matrices)
449  typename MatrixType::InnerIterator it_spmat(spmat, j);
450  for(; it_spmat; ++it_spmat)
451  {
452  StorageIndex idx = 0; // Position of this block in the column block
453  StorageIndex bi = innerToBlock(it_spmat.index()); // Index of the current nonzero block
454  // Go to the inner block where this element belongs to
455  while(bi > m_indices[m_outerIndex[bj]+idx]) ++idx; // Not expensive for ordered blocks
456  StorageIndex idxVal;// Get the right position in the array of values for this element
457  if(m_blockSize == Dynamic)
458  {
459  // Offset from all blocks before ...
460  idxVal = m_blockPtr[m_outerIndex[bj]+idx];
461  // ... and offset inside the block
462  idxVal += (j - blockOuterIndex(bj)) * blockOuterSize(bj) + it_spmat.index() - m_innerOffset[bi];
463  }
464  else
465  {
466  // All blocks before
467  idxVal = (m_outerIndex[bj] + idx) * m_blockSize * m_blockSize;
468  // inside the block
469  idxVal += (j - blockOuterIndex(bj)) * m_blockSize + (it_spmat.index()%m_blockSize);
470  }
471  // Insert the value
472  m_values[idxVal] = it_spmat.value();
473  } // end of this column
474  } // end of this block
475  } // end of this outer block
476 
477  return *this;
478  }
479 
497  template<typename MatrixType>
498  void setBlockStructure(const MatrixType& blockPattern)
499  {
500  resize(blockPattern.rows(), blockPattern.cols());
501  reserve(blockPattern.nonZeros());
502 
503  // Browse the block pattern and set up the various pointers
504  m_outerIndex[0] = 0;
505  if(m_blockSize == Dynamic) m_blockPtr[0] = 0;
506  for(StorageIndex nz = 0; nz < m_nonzeros; ++nz) m_values[nz] = Scalar(0);
507  for(StorageIndex bj = 0; bj < m_outerBSize; ++bj)
508  {
509  //Browse each outer block
510 
511  //First, copy and save the indices of nonzero blocks
512  //FIXME : find a way to avoid this ...
513  std::vector<int> nzBlockIdx;
514  typename MatrixType::InnerIterator it(blockPattern, bj);
515  for(; it; ++it)
516  {
517  nzBlockIdx.push_back(it.index());
518  }
519  std::sort(nzBlockIdx.begin(), nzBlockIdx.end());
520 
521  // Now, fill block indices and (eventually) pointers to blocks
522  for(StorageIndex idx = 0; idx < nzBlockIdx.size(); ++idx)
523  {
524  StorageIndex offset = m_outerIndex[bj]+idx; // offset in m_indices
525  m_indices[offset] = nzBlockIdx[idx];
526  if(m_blockSize == Dynamic)
527  m_blockPtr[offset] = m_blockPtr[offset-1] + blockInnerSize(nzBlockIdx[idx]) * blockOuterSize(bj);
528  // There is no blockPtr for fixed-size blocks... not needed !???
529  }
530  // Save the pointer to the next outer block
531  m_outerIndex[bj+1] = m_outerIndex[bj] + nzBlockIdx.size();
532  }
533  }
534 
538  inline void resize(Index brow, Index bcol)
539  {
540  m_innerBSize = IsColMajor ? brow : bcol;
541  m_outerBSize = IsColMajor ? bcol : brow;
542  }
543 
549  inline void setBlockSize(Index blockSize)
550  {
551  m_blockSize = blockSize;
552  }
553 
563  inline void setBlockLayout(const VectorXi& rowBlocks, const VectorXi& colBlocks)
564  {
565  const VectorXi& innerBlocks = IsColMajor ? rowBlocks : colBlocks;
566  const VectorXi& outerBlocks = IsColMajor ? colBlocks : rowBlocks;
567  eigen_assert(m_innerBSize == innerBlocks.size() && "CHECK THE NUMBER OF ROW OR COLUMN BLOCKS");
568  eigen_assert(m_outerBSize == outerBlocks.size() && "CHECK THE NUMBER OF ROW OR COLUMN BLOCKS");
569  m_outerBSize = outerBlocks.size();
570  // starting index of blocks... cumulative sums
571  m_innerOffset = new StorageIndex[m_innerBSize+1];
572  m_outerOffset = new StorageIndex[m_outerBSize+1];
573  m_innerOffset[0] = 0;
574  m_outerOffset[0] = 0;
575  std::partial_sum(&innerBlocks[0], &innerBlocks[m_innerBSize-1]+1, &m_innerOffset[1]);
576  std::partial_sum(&outerBlocks[0], &outerBlocks[m_outerBSize-1]+1, &m_outerOffset[1]);
577 
578  // Compute the total number of nonzeros
579  m_nonzeros = 0;
580  for(StorageIndex bj = 0; bj < m_outerBSize; ++bj)
581  for(StorageIndex bi = 0; bi < m_innerBSize; ++bi)
582  m_nonzeros += outerBlocks[bj] * innerBlocks[bi];
583 
584  }
585 
596  inline void reserve(const Index nonzerosblocks)
597  {
598  eigen_assert((m_innerBSize != 0 && m_outerBSize != 0) &&
599  "TRYING TO RESERVE ZERO-SIZE MATRICES, CALL resize() first");
600 
601  //FIXME Should free if already allocated
602  m_outerIndex = new StorageIndex[m_outerBSize+1];
603 
604  m_nonzerosblocks = nonzerosblocks;
605  if(m_blockSize != Dynamic)
606  {
607  m_nonzeros = nonzerosblocks * (m_blockSize * m_blockSize);
608  m_blockPtr = 0;
609  }
610  else
611  {
612  // m_nonzeros is already computed in setBlockLayout()
613  m_blockPtr = new StorageIndex[m_nonzerosblocks+1];
614  }
615  m_indices = new StorageIndex[m_nonzerosblocks+1];
616  m_values = new Scalar[m_nonzeros];
617  }
618 
619 
630  template<typename InputIterator>
631  void setFromTriplets(const InputIterator& begin, const InputIterator& end)
632  {
633  eigen_assert((m_innerBSize!=0 && m_outerBSize !=0) && "ZERO BLOCKS, PLEASE CALL resize() before");
634 
635  /* First, sort the triplet list
636  * FIXME This can be unnecessarily expensive since only the inner indices have to be sorted
637  * The best approach is like in SparseMatrix::setFromTriplets()
638  */
640  std::sort(begin, end, tripletcomp);
641 
642  /* Count the number of rows and column blocks,
643  * and the number of nonzero blocks per outer dimension
644  */
645  VectorXi rowBlocks(m_innerBSize); // Size of each block row
646  VectorXi colBlocks(m_outerBSize); // Size of each block column
647  rowBlocks.setZero(); colBlocks.setZero();
648  VectorXi nzblock_outer(m_outerBSize); // Number of nz blocks per outer vector
649  VectorXi nz_outer(m_outerBSize); // Number of nz per outer vector...for variable-size blocks
650  nzblock_outer.setZero();
651  nz_outer.setZero();
652  for(InputIterator it(begin); it !=end; ++it)
653  {
654  eigen_assert(it->row() >= 0 && it->row() < this->blockRows() && it->col() >= 0 && it->col() < this->blockCols());
655  eigen_assert((it->value().rows() == it->value().cols() && (it->value().rows() == m_blockSize))
656  || (m_blockSize == Dynamic));
657 
658  if(m_blockSize == Dynamic)
659  {
660  eigen_assert((rowBlocks[it->row()] == 0 || rowBlocks[it->row()] == it->value().rows()) &&
661  "NON CORRESPONDING SIZES FOR ROW BLOCKS");
662  eigen_assert((colBlocks[it->col()] == 0 || colBlocks[it->col()] == it->value().cols()) &&
663  "NON CORRESPONDING SIZES FOR COLUMN BLOCKS");
664  rowBlocks[it->row()] =it->value().rows();
665  colBlocks[it->col()] = it->value().cols();
666  }
667  nz_outer(IsColMajor ? it->col() : it->row()) += it->value().rows() * it->value().cols();
668  nzblock_outer(IsColMajor ? it->col() : it->row())++;
669  }
670  // Allocate member arrays
671  if(m_blockSize == Dynamic) setBlockLayout(rowBlocks, colBlocks);
672  StorageIndex nzblocks = nzblock_outer.sum();
673  reserve(nzblocks);
674 
675  // Temporary markers
676  VectorXi block_id(m_outerBSize); // To be used as a block marker during insertion
677 
678  // Setup outer index pointers and markers
679  m_outerIndex[0] = 0;
680  if (m_blockSize == Dynamic) m_blockPtr[0] = 0;
681  for(StorageIndex bj = 0; bj < m_outerBSize; ++bj)
682  {
683  m_outerIndex[bj+1] = m_outerIndex[bj] + nzblock_outer(bj);
684  block_id(bj) = m_outerIndex[bj];
685  if(m_blockSize==Dynamic)
686  {
687  m_blockPtr[m_outerIndex[bj+1]] = m_blockPtr[m_outerIndex[bj]] + nz_outer(bj);
688  }
689  }
690 
691  // Fill the matrix
692  for(InputIterator it(begin); it!=end; ++it)
693  {
694  StorageIndex outer = IsColMajor ? it->col() : it->row();
695  StorageIndex inner = IsColMajor ? it->row() : it->col();
696  m_indices[block_id(outer)] = inner;
697  StorageIndex block_size = it->value().rows()*it->value().cols();
698  StorageIndex nz_marker = blockPtr(block_id[outer]);
699  memcpy(&(m_values[nz_marker]), it->value().data(), block_size * sizeof(Scalar));
700  if(m_blockSize == Dynamic)
701  {
702  m_blockPtr[block_id(outer)+1] = m_blockPtr[block_id(outer)] + block_size;
703  }
704  block_id(outer)++;
705  }
706 
707  // An alternative when the outer indices are sorted...no need to use an array of markers
708 // for(Index bcol = 0; bcol < m_outerBSize; ++bcol)
709 // {
710 // Index id = 0, id_nz = 0, id_nzblock = 0;
711 // for(InputIterator it(begin); it!=end; ++it)
712 // {
713 // while (id<bcol) // one pass should do the job unless there are empty columns
714 // {
715 // id++;
716 // m_outerIndex[id+1]=m_outerIndex[id];
717 // }
718 // m_outerIndex[id+1] += 1;
719 // m_indices[id_nzblock]=brow;
720 // Index block_size = it->value().rows()*it->value().cols();
721 // m_blockPtr[id_nzblock+1] = m_blockPtr[id_nzblock] + block_size;
722 // id_nzblock++;
723 // memcpy(&(m_values[id_nz]),it->value().data(), block_size*sizeof(Scalar));
724 // id_nz += block_size;
725 // }
726 // while(id < m_outerBSize-1) // Empty columns at the end
727 // {
728 // id++;
729 // m_outerIndex[id+1]=m_outerIndex[id];
730 // }
731 // }
732  }
733 
734 
738  inline Index rows() const
739  {
740 // return blockRows();
741  return (IsColMajor ? innerSize() : outerSize());
742  }
743 
747  inline Index cols() const
748  {
749 // return blockCols();
750  return (IsColMajor ? outerSize() : innerSize());
751  }
752 
753  inline Index innerSize() const
754  {
755  if(m_blockSize == Dynamic) return m_innerOffset[m_innerBSize];
756  else return (m_innerBSize * m_blockSize) ;
757  }
758 
759  inline Index outerSize() const
760  {
761  if(m_blockSize == Dynamic) return m_outerOffset[m_outerBSize];
762  else return (m_outerBSize * m_blockSize) ;
763  }
765  inline Index blockRows() const
766  {
767  return (IsColMajor ? m_innerBSize : m_outerBSize);
768  }
770  inline Index blockCols() const
771  {
772  return (IsColMajor ? m_outerBSize : m_innerBSize);
773  }
774 
775  inline Index outerBlocks() const { return m_outerBSize; }
776  inline Index innerBlocks() const { return m_innerBSize; }
777 
779  inline Index outerToBlock(Index outer) const
780  {
781  eigen_assert(outer < outerSize() && "OUTER INDEX OUT OF BOUNDS");
782 
783  if(m_blockSize != Dynamic)
784  return (outer / m_blockSize); // Integer division
785 
786  StorageIndex b_outer = 0;
787  while(m_outerOffset[b_outer] <= outer) ++b_outer;
788  return b_outer - 1;
789  }
791  inline Index innerToBlock(Index inner) const
792  {
793  eigen_assert(inner < innerSize() && "OUTER INDEX OUT OF BOUNDS");
794 
795  if(m_blockSize != Dynamic)
796  return (inner / m_blockSize); // Integer division
797 
798  StorageIndex b_inner = 0;
799  while(m_innerOffset[b_inner] <= inner) ++b_inner;
800  return b_inner - 1;
801  }
802 
807  {
808  eigen_assert(brow < blockRows() && "BLOCK ROW INDEX OUT OF BOUNDS");
809  eigen_assert(bcol < blockCols() && "BLOCK nzblocksFlagCOLUMN OUT OF BOUNDS");
810 
811  StorageIndex rsize = IsColMajor ? blockInnerSize(brow): blockOuterSize(bcol);
812  StorageIndex csize = IsColMajor ? blockOuterSize(bcol) : blockInnerSize(brow);
813  StorageIndex inner = IsColMajor ? brow : bcol;
814  StorageIndex outer = IsColMajor ? bcol : brow;
815  StorageIndex offset = m_outerIndex[outer];
816  while(offset < m_outerIndex[outer+1] && m_indices[offset] != inner)
817  offset++;
818  if(m_indices[offset] == inner)
819  {
820  return Map<BlockScalar>(&(m_values[blockPtr(offset)]), rsize, csize);
821  }
822  else
823  {
824  //FIXME the block does not exist, Insert it !!!!!!!!!
825  eigen_assert("DYNAMIC INSERTION IS NOT YET SUPPORTED");
826  }
827  }
828 
833  {
834  eigen_assert(brow < blockRows() && "BLOCK ROW INDEX OUT OF BOUNDS");
835  eigen_assert(bcol < blockCols() && "BLOCK COLUMN OUT OF BOUNDS");
836 
837  StorageIndex rsize = IsColMajor ? blockInnerSize(brow): blockOuterSize(bcol);
838  StorageIndex csize = IsColMajor ? blockOuterSize(bcol) : blockInnerSize(brow);
839  StorageIndex inner = IsColMajor ? brow : bcol;
840  StorageIndex outer = IsColMajor ? bcol : brow;
841  StorageIndex offset = m_outerIndex[outer];
842  while(offset < m_outerIndex[outer+1] && m_indices[offset] != inner) offset++;
843  if(m_indices[offset] == inner)
844  {
845  return Map<const BlockScalar> (&(m_values[blockPtr(offset)]), rsize, csize);
846  }
847  else
848 // return BlockScalar::Zero(rsize, csize);
849  eigen_assert("NOT YET SUPPORTED");
850  }
851 
852  // Block Matrix times vector product
853  template<typename VecType>
855  {
857  }
858 
860  inline Index nonZerosBlocks() const { return m_nonzerosblocks; }
862  inline Index nonZeros() const { return m_nonzeros; }
863 
864  inline BlockScalarReturnType *valuePtr() {return static_cast<BlockScalarReturnType *>(m_values);}
865 // inline Scalar *valuePtr(){ return m_values; }
866  inline StorageIndex *innerIndexPtr() {return m_indices; }
867  inline const StorageIndex *innerIndexPtr() const {return m_indices; }
868  inline StorageIndex *outerIndexPtr() {return m_outerIndex; }
869  inline const StorageIndex* outerIndexPtr() const {return m_outerIndex; }
870 
872  inline bool isCompressed() const {return true;}
876  inline Index blockRowsIndex(Index bi) const
877  {
878  return IsColMajor ? blockInnerIndex(bi) : blockOuterIndex(bi);
879  }
880 
884  inline Index blockColsIndex(Index bj) const
885  {
886  return IsColMajor ? blockOuterIndex(bj) : blockInnerIndex(bj);
887  }
888 
889  inline Index blockOuterIndex(Index bj) const
890  {
891  return (m_blockSize == Dynamic) ? m_outerOffset[bj] : (bj * m_blockSize);
892  }
893  inline Index blockInnerIndex(Index bi) const
894  {
895  return (m_blockSize == Dynamic) ? m_innerOffset[bi] : (bi * m_blockSize);
896  }
897 
898  // Not needed ???
899  inline Index blockInnerSize(Index bi) const
900  {
901  return (m_blockSize == Dynamic) ? (m_innerOffset[bi+1] - m_innerOffset[bi]) : m_blockSize;
902  }
903  inline Index blockOuterSize(Index bj) const
904  {
905  return (m_blockSize == Dynamic) ? (m_outerOffset[bj+1]- m_outerOffset[bj]) : m_blockSize;
906  }
907 
911  class InnerIterator; // Browse column by column
912 
916  class BlockInnerIterator; // Browse block by block
917 
918  friend std::ostream & operator << (std::ostream & s, const BlockSparseMatrix& m)
919  {
920  for (StorageIndex j = 0; j < m.outerBlocks(); ++j)
921  {
922  BlockInnerIterator itb(m, j);
923  for(; itb; ++itb)
924  {
925  s << "("<<itb.row() << ", " << itb.col() << ")\n";
926  s << itb.value() <<"\n";
927  }
928  }
929  s << std::endl;
930  return s;
931  }
932 
936  Index blockPtr(Index id) const
937  {
938  if(m_blockSize == Dynamic) return m_blockPtr[id];
939  else return id * m_blockSize * m_blockSize;
940  //return blockDynIdx(id, typename internal::conditional<(BlockSize==Dynamic), internal::true_type, internal::false_type>::type());
941  }
942 
943 
944  protected:
945 // inline Index blockDynIdx(Index id, internal::true_type) const
946 // {
947 // return m_blockPtr[id];
948 // }
949 // inline Index blockDynIdx(Index id, internal::false_type) const
950 // {
951 // return id * BlockSize * BlockSize;
952 // }
953 
954  // To be implemented
955  // Insert a block at a particular location... need to make a room for that
956  Map<BlockScalar> insert(Index brow, Index bcol);
957 
958  Index m_innerBSize; // Number of block rows
959  Index m_outerBSize; // Number of block columns
960  StorageIndex *m_innerOffset; // Starting index of each inner block (size m_innerBSize+1)
961  StorageIndex *m_outerOffset; // Starting index of each outer block (size m_outerBSize+1)
962  Index m_nonzerosblocks; // Total nonzeros blocks (lower than m_innerBSize x m_outerBSize)
963  Index m_nonzeros; // Total nonzeros elements
964  Scalar *m_values; //Values stored block column after block column (size m_nonzeros)
965  StorageIndex *m_blockPtr; // Pointer to the beginning of each block in m_values, size m_nonzeroblocks ... null for fixed-size blocks
966  StorageIndex *m_indices; //Inner block indices, size m_nonzerosblocks ... OK
967  StorageIndex *m_outerIndex; // Starting pointer of each block column in m_indices (size m_outerBSize)... OK
968  Index m_blockSize; // Size of a block for fixed-size blocks, otherwise -1
969 };
970 
971 template<typename _Scalar, int _BlockAtCompileTime, int _Options, typename _StorageIndex>
972 class BlockSparseMatrix<_Scalar, _BlockAtCompileTime, _Options, _StorageIndex>::BlockInnerIterator
973 {
974  public:
975 
976  enum{
977  Flags = _Options
978  };
979 
981  : m_mat(mat),m_outer(outer),
982  m_id(mat.m_outerIndex[outer]),
983  m_end(mat.m_outerIndex[outer+1])
984  {
985  }
986 
987  inline BlockInnerIterator& operator++() {m_id++; return *this; }
988 
989  inline const Map<const BlockScalar> value() const
990  {
991  return Map<const BlockScalar>(&(m_mat.m_values[m_mat.blockPtr(m_id)]),
992  rows(),cols());
993  }
995  {
996  return Map<BlockScalar>(&(m_mat.m_values[m_mat.blockPtr(m_id)]),
997  rows(),cols());
998  }
999  // Block inner index
1000  inline Index index() const {return m_mat.m_indices[m_id]; }
1001  inline Index outer() const { return m_outer; }
1002  // block row index
1003  inline Index row() const {return index(); }
1004  // block column index
1005  inline Index col() const {return outer(); }
1006  // FIXME Number of rows in the current block
1007  inline Index rows() const { return (m_mat.m_blockSize==Dynamic) ? (m_mat.m_innerOffset[index()+1] - m_mat.m_innerOffset[index()]) : m_mat.m_blockSize; }
1008  // Number of columns in the current block ...
1009  inline Index cols() const { return (m_mat.m_blockSize==Dynamic) ? (m_mat.m_outerOffset[m_outer+1]-m_mat.m_outerOffset[m_outer]) : m_mat.m_blockSize;}
1010  inline operator bool() const { return (m_id < m_end); }
1011 
1012  protected:
1017 };
1018 
1019 template<typename _Scalar, int _BlockAtCompileTime, int _Options, typename _StorageIndex>
1020 class BlockSparseMatrix<_Scalar, _BlockAtCompileTime, _Options, _StorageIndex>::InnerIterator
1021 {
1022  public:
1024  : m_mat(mat),m_outerB(mat.outerToBlock(outer)),m_outer(outer),
1025  itb(mat, mat.outerToBlock(outer)),
1026  m_offset(outer - mat.blockOuterIndex(m_outerB))
1027  {
1028  if (itb)
1029  {
1030  m_id = m_mat.blockInnerIndex(itb.index());
1031  m_start = m_id;
1032  m_end = m_mat.blockInnerIndex(itb.index()+1);
1033  }
1034  }
1036  {
1037  m_id++;
1038  if (m_id >= m_end)
1039  {
1040  ++itb;
1041  if (itb)
1042  {
1043  m_id = m_mat.blockInnerIndex(itb.index());
1044  m_start = m_id;
1045  m_end = m_mat.blockInnerIndex(itb.index()+1);
1046  }
1047  }
1048  return *this;
1049  }
1050  inline const Scalar& value() const
1051  {
1052  return itb.value().coeff(m_id - m_start, m_offset);
1053  }
1054  inline Scalar& valueRef()
1055  {
1056  return itb.valueRef().coeff(m_id - m_start, m_offset);
1057  }
1058  inline Index index() const { return m_id; }
1059  inline Index outer() const {return m_outer; }
1060  inline Index col() const {return outer(); }
1061  inline Index row() const { return index();}
1062  inline operator bool() const
1063  {
1064  return itb;
1065  }
1066  protected:
1070  BlockInnerIterator itb; // Iterator through the blocks
1071  const Index m_offset; // Position of this column in the block
1072  Index m_start; // starting inner index of this block
1073  Index m_id; // current inner index in the block
1074  Index m_end; // starting inner index of the next block
1075 
1076 };
1077 } // end namespace Eigen
1078 
1079 #endif // EIGEN_SPARSEBLOCKMATRIX_H
Matrix3f m
void sparse_time_dense_product(const SparseLhsType &lhs, const DenseRhsType &rhs, DenseResType &res, const AlphaType &alpha)
Scalar coeff(Index bi) const
void setBlockSize(Index blockSize)
set the block size at runtime for fixed-size block layout
void setFromTriplets(const InputIterator &begin, const InputIterator &end)
Fill values in a matrix from a triplet list.
BlockSparseTimeDenseProduct< BlockSparseMatrix, VecType > operator*(const VecType &lhs) const
SCALAR Scalar
Definition: bench_gemm.cpp:46
Index innerToBlock(Index inner) const
Matrix< RealScalar, _BlockAtCompileTime, _BlockAtCompileTime, IsColMajor ? ColMajor :RowMajor > BlockRealScalar
const BlockSparseMatrixT & m_spblockmat
A insert(1, 2)=0
Index blockOuterSize(Index bj) const
Index blockRowsIndex(Index bi) const
Scalar * b
Definition: benchVecAdd.cpp:17
void setBlockStructure(const MatrixType &blockPattern)
Set the nonzero block pattern of the matrix.
A versatible sparse matrix representation.
Definition: SparseMatrix.h:96
const BlockSparseMatrixT & m_spblockmat
BlockSparseMatrix< Scalar, BlockSize, IsColMajor ? ColMajor :RowMajor, StorageIndex > PlainObject
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
const int InnerRandomAccessPattern
Definition: SparseUtil.h:48
BlockSparseMatrixView(const BlockSparseMatrixT &spblockmat)
Index blockInnerIndex(Index bi) const
const BlockSparseMatrix< _Scalar, _BlockAtCompileTime, _Options, StorageIndex > & m_mat
const unsigned int LvalueBit
Definition: Constants.h:144
InnerIterator(const BlockSparseMatrix &mat, Index outer)
BlockSparseMatrix(const MatrixType &spmat)
Constructor from a sparse matrix.
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy y set format x g set format y g set format x2 g set format y2 g set format z g set angles radians set nogrid set key title set key left top Right noreverse box linetype linewidth samplen spacing width set nolabel set noarrow set nologscale set logscale x set set pointsize set encoding default set nopolar set noparametric set set set set surface set nocontour set clabel set mapping cartesian set nohidden3d set cntrparam order set cntrparam linear set cntrparam levels auto set cntrparam points set size set set xzeroaxis lt lw set x2zeroaxis lt lw set yzeroaxis lt lw set y2zeroaxis lt lw set tics in set ticslevel set tics set mxtics default set mytics default set mx2tics default set my2tics default set xtics border mirror norotate autofreq set ytics border mirror norotate autofreq set ztics border nomirror norotate autofreq set nox2tics set noy2tics set timestamp bottom norotate offset
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
StorageIndex * outerIndexPtr()
BlockScalarReturnType * valuePtr()
const StorageIndex * outerIndexPtr() const
MatrixXf MatrixType
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
const Map< const BlockScalar > value() const
#define EIGEN_STATIC_ASSERT(CONDITION, MSG)
Definition: StaticAssert.h:127
BlockSparseMatrix(Index brow, Index bcol)
Construct and resize.
EIGEN_CONSTEXPR Index first(const T &x) EIGEN_NOEXCEPT
BlockInnerIterator(const BlockSparseMatrix &mat, const Index outer)
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:39
const unsigned int RowMajorBit
Definition: Constants.h:66
Scalar coeff(Index row, Index col)
static const Similarity3 id
const BlockSparseMatrixT & m_spblockmat
NumTraits< Scalar >::Real RealScalar
StorageIndex * innerIndexPtr()
internal::conditional< _BlockAtCompileTime==Dynamic, Scalar, BlockScalar >::type BlockScalarReturnType
void reserve(const Index nonzerosblocks)
Allocate the internal array of pointers to blocks and their inner indices.
Ref< typename BlockSparseMatrixT::BlockRealScalar > RealScalar
BlockSparseMatrix & operator=(BlockSparseMatrix other)
Scalar coeff(Index bi, Index j) const
Base class of any sparse matrices or sparse expressions.
BlockSparseMatrix & operator=(const MatrixType &spmat)
Assignment from a sparse matrix with the same storage order.
const VectorType & m_vec
bool isCompressed() const
for compatibility purposes with the SparseMatrix class
void setBlockLayout(const VectorXi &rowBlocks, const VectorXi &colBlocks)
Set the row and column block layouts,.
m row(1)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
BlockSparseMatrixT::Index Index
#define eigen_assert(x)
Definition: Macros.h:1037
friend void swap(BlockSparseMatrix &first, BlockSparseMatrix &second)
BlockSparseMatrix(const BlockSparseMatrix &other)
Copy-constructor.
RealScalar alpha
RealScalar s
Ref< typename BlockSparseMatrixT::BlockScalar > Scalar
Index outerToBlock(Index outer) const
Index blockOuterIndex(Index bj) const
Scalar coeffRef(Index bi, Index j)
Matrix< Scalar, _BlockAtCompileTime, _BlockAtCompileTime, IsColMajor ? ColMajor :RowMajor > BlockScalar
Ref< const Matrix< typename BlockSparseMatrixT::Scalar,(RowsAtCompileTime==1)? 1 :BlockSize,(ColsAtCompileTime==1)? 1 :BlockSize > > Scalar
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
Ref< Matrix< typename BlockSparseMatrixT::RealScalar, BlockSparseMatrixT::BlockSize, BlockSparseMatrixT::BlockSize > > RealScalar
Iterator::value_type Triplet
BlockSparseMatrixT::Index Index
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:281
Map< const BlockScalar > coeff(Index brow, Index bcol) const
const unsigned int NestByRefBit
Definition: Constants.h:169
BlockSparseMatrixT::Index Index
Index blockPtr(Index id) const
BlockVectorReturn(const BlockSparseMatrixT &spblockmat, VectorType &vec)
static EIGEN_DEPRECATED const end_t end
void scaleAndAddTo(Dest &dest, const typename Rhs::Scalar &alpha) const
Ref< Matrix< typename BlockSparseMatrixT::Scalar, BlockSparseMatrixT::BlockSize, BlockSparseMatrixT::BlockSize > > Scalar
A versatile sparse matrix representation where each element is a block.
m col(1)
const StorageIndex * innerIndexPtr() const
const int Dynamic
Definition: Constants.h:22
Index blockInnerSize(Index bi) const
Ref< Matrix< typename VectorType::Scalar, RowsAtCompileTime, ColsAtCompileTime > > Scalar
Generic expression where a coefficient-wise unary operator is applied to an expression.
Definition: CwiseUnaryOp.h:55
The matrix class, also used for vectors and row-vectors.
Index blockColsIndex(Index bj) const
v resize(3)
BlockVectorView(const BlockSparseMatrixT &spblockmat, const VectorType &vec)
std::ostream & operator<<(std::ostream &s, const Packet16c &v)
std::ptrdiff_t j
Scalar coeffRef(Index row, Index col)
bool operator()(const Triplet &a, const Triplet &b)
internal::ref_selector< BlockSparseMatrix< _Scalar, _BlockAtCompileTime, _Options, _StorageIndex > >::type Nested
void resize(Index brow, Index bcol)
Set the number of rows and columns blocks.
Definition: pytypes.h:1370
InnerIterator(const BlockSparseMatrixView &mat, Index outer)
void swap(scoped_array< T > &a, scoped_array< T > &b)
Definition: Memory.h:709
An InnerIterator allows to loop over the element of any matrix expression.
Definition: CoreIterators.h:33
Ref< BlockScalar > coeffRef(Index brow, Index bcol)


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:33:59