So-Bogus
A c++ sparse block matrix library aimed at Second Order cone problems
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Pages
SparseBlockMatrixBase.hpp
1 /*
2  * This file is part of bogus, a C++ sparse block matrix library.
3  *
4  * Copyright 2013 Gilles Daviet <gdaviet@gmail.com>
5  *
6  * This Source Code Form is subject to the terms of the Mozilla Public
7  * License, v. 2.0. If a copy of the MPL was not distributed with this
8  * file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 */
10 
11 
12 #ifndef BOGUS_SPARSEBLOCKMATRIX_BASE_HPP
13 #define BOGUS_SPARSEBLOCKMATRIX_BASE_HPP
14 
15 #include "BlockMatrixBase.hpp"
16 #include "Expressions.hpp"
17 #include "DynamicExpressions.hpp"
18 
19 #include "SparseBlockIndex.hpp"
20 #include "CompressedSparseBlockIndex.hpp"
21 
22 #include "../Utils/CppTools.hpp"
23 #include "../Utils/Lock.hpp"
24 
25 namespace bogus
26 {
27 
28 template < typename BlockT, int Flags > class SparseBlockMatrix ;
29 template < bool Symmetric > struct SparseBlockMatrixFinalizer ;
30 template < typename Derived, bool Major > struct SparseBlockIndexGetter ;
31 
33 
35 template < typename Derived >
36 class SparseBlockMatrixBase : public BlockMatrixBase< Derived >
37 {
38 
39 public:
40  // Convenient typedefs and using directives
41 
44 
45  typedef typename Traits::BlockPtr BlockPtr ;
46  typedef typename Base::Index Index ;
47 
48 
49  typedef typename Traits::MajorIndexType MajorIndexType ;
50 
51  typedef typename Base::BlockType BlockType ;
52  typedef typename Base::Scalar Scalar ;
53 
54  typedef typename MajorIndexType::InnerIterator InnerIterator ;
55 
58 
59  // Minor index is always uncompressed, as the blocks cannot be contiguous
60  // For a symmetric matrix, it does not store diagonal block in the minor and transpose index
62  // Transpose index is compressed for perf, as we always create it in a compressed-compatible way
64 
65  typedef typename TypeSwapIf< Traits::is_col_major, MajorIndexType, MinorIndexType >::First RowIndexType ;
66  typedef typename TypeSwapIf< Traits::is_col_major, MajorIndexType, MinorIndexType >::Second ColIndexType ;
67 
68  // Canonical type for a mutable matrix with different block type
69  template < typename OtherBlockType, bool PreserveSymmetry = true, bool SwitchDirection = false >
70  struct MutableImpl
71  {
72  typedef SparseBlockMatrix< OtherBlockType,
73  Traits::flags & ~flags::UNCOMPRESSED > Type ;
74  } ;
75  template < typename OtherBlockType >
76  struct MutableImpl< OtherBlockType, false, false >
77  {
78  typedef SparseBlockMatrix< OtherBlockType,
79  Traits::flags & ~flags::UNCOMPRESSED & ~flags::SYMMETRIC > Type ;
80  } ;
81  template < typename OtherBlockType >
82  struct MutableImpl< OtherBlockType, true, true >
83  {
84  typedef SparseBlockMatrix< OtherBlockType,
85  (Traits::flags & ~flags::UNCOMPRESSED &
86  ~flags::COL_MAJOR) | ((~Traits::flags)&flags::COL_MAJOR) > Type ;
87  } ;
88  template < typename OtherBlockType >
89  struct MutableImpl< OtherBlockType, false, true >
90  {
91  typedef SparseBlockMatrix< OtherBlockType,
92  (Traits::flags & ~flags::UNCOMPRESSED & ~flags::SYMMETRIC&
93  ~flags::COL_MAJOR) | ((~Traits::flags)&flags::COL_MAJOR) > Type ;
94  } ;
96 
97  typedef typename Base::ConstTransposeReturnType ConstTransposeReturnType ;
98 
99  typedef typename BlockTraits< BlockType >::TransposeStorageType TransposeBlockType ;
100  typedef typename ResizableSequenceContainer< TransposeBlockType >::Type TransposeArrayType ;
101 
102  using Base::rows ;
103  using Base::cols ;
104  using Base::blocks ;
105  using Base::derived ;
106  using Base::block ;
107  using Base::diagonal ;
108  using Base::InvalidBlockPtr ;
109 
110  // Public API
111 
114 
116 
119  void setRows( const Index nBlocks, const unsigned* rowsPerBlock ) ;
121  void setRows( const std::vector< unsigned > &rowsPerBlock )
122  {
123  setRows( rowsPerBlock.size(), &rowsPerBlock[0] ) ;
124  }
126  void setRows( const Index nBlocks, const Index rowsPerBlock )
127  {
128  setRows( std::vector< unsigned >( nBlocks, rowsPerBlock ) ) ;
129  }
131  void setRows( const Index nBlocks )
132  {
133  assert( Base::has_fixed_rows_blocks ) ;
134  setRows( nBlocks, BlockType::RowsAtCompileTime ) ;
135  }
136 
138 
141  void setCols( const Index nBlocks, const unsigned* colsPerBlock ) ;
143  void setCols( const std::vector< unsigned > &colsPerBlock )
144  {
145  setCols( colsPerBlock.size(), &colsPerBlock[0] ) ;
146  }
148  void setCols( const Index nBlocks, const Index colsPerBlock )
149  {
150  setCols( std::vector< unsigned >( nBlocks, colsPerBlock ) ) ;
151  }
153  void setCols( const Index nBlocks )
154  {
155  assert( Base::has_fixed_cols_blocks ) ;
156  setCols( nBlocks, BlockType::ColsAtCompileTime ) ;
157  }
158 
159  Index rowsOfBlocks() const { return rowMajorIndex().outerSize() ; }
160  Index colsOfBlocks() const { return colMajorIndex().outerSize() ; }
161 
162  Index blockRows( Index row ) const { return rowOffsets()[ row + 1 ] - rowOffsets()[ row ] ; }
163  Index blockCols( Index col ) const { return colOffsets()[ col + 1 ] - colOffsets()[ col ] ; }
164 
166 
169 
171  void reserve( std::size_t nBlocks )
172  {
173  m_blocks.reserve( nBlocks ) ;
174  m_majorIndex.reserve( nBlocks ) ;
175  }
176 
178 
183  BlockType& insertBack( Index row, Index col )
184  {
185  if( Traits::is_col_major )
186  return insertByOuterInner< true >( col, row ) ;
187  else
188  return insertByOuterInner< true >( row, col ) ;
189  }
190 
192  BlockType& insertBackAndResize( Index row, Index col ) ;
193 
195 
208  BlockType& insert( Index row, Index col )
209  {
210  BOGUS_STATIC_ASSERT( !Traits::is_compressed, UNORDERED_INSERTION_WITH_COMPRESSED_INDEX ) ;
211  if( Traits::is_col_major )
212  return insertByOuterInner< false >( col, row ) ;
213  else
214  return insertByOuterInner< false >( row, col ) ;
215  }
216 
218  BlockType& insertAndResize( Index row, Index col ) ;
219 
221 
223  template< bool Ordered >
224  BlockType& insertByOuterInner( Index outer, Index inner ) ;
225 
228  void finalize() ;
229 
231  void clear() ;
233  Derived& setZero() { clear() ; return derived() ; }
235  Derived& setIdentity() ;
237  Derived& setBlocksToZero() ;
238 
240  std::size_t nBlocks() const { return blocks().size() ; }
241  std::size_t size() const { return nBlocks() ; }
242 
244  bool empty() const { return 0 == nBlocks() ; }
245 
246  const BlockType& block( BlockPtr ptr ) const
247  {
248  return m_blocks[ ptr ] ;
249  }
250 
251  BlockType& block( BlockPtr ptr )
252  {
253  return m_blocks[ ptr ] ;
254  }
255 
257  BlockPtr blockPtr( Index row, Index col ) const ;
259  BlockPtr diagonalBlockPtr( Index row ) const ;
260 
262  template <bool ColWise, typename Func>
263  void eachBlockOf( const Index outer, Func func ) const ;
264 
266 
269 
271 
275  bool computeMinorIndex() ;
276 
278 
281  void cacheTranspose() ;
283  bool transposeCached() const { return m_transposeIndex.valid ; }
284 
286 
295  InnerIterator innerIterator( Index outer ) const
296  {
297  return InnerIterator( m_majorIndex, outer ) ;
298  }
299 
300  const MajorIndexType& majorIndex() const
301  {
302  return m_majorIndex ;
303  }
304  const MinorIndexType& minorIndex() const
305  {
306  return m_minorIndex ;
307  }
308  const TransposeIndexType& transposeIndex() const
309  {
310  return m_transposeIndex ;
311  }
312  const ColIndexType& colMajorIndex() const ;
313  const RowIndexType& rowMajorIndex() const ;
314 
316 
319 
321  template < bool Transpose, typename OtherDerived >
322  Derived& assign ( const SparseBlockMatrixBase< OtherDerived > &source, const Scalar scale = 1 ) ;
323 
324  template < typename OtherDerived >
325  Derived& operator= ( const BlockObjectBase< OtherDerived > &source )
326  {
327  Evaluator< OtherDerived > rhs ( source.derived() ) ;
328  return assign< OtherDerived::is_transposed >( *rhs ) ;
329  }
330 
331  Derived& operator= ( const SparseBlockMatrixBase &source ) ;
332 
333  template < typename LhsT, typename RhsT >
334  Derived& operator= ( const Product< LhsT, RhsT > &prod ) ;
335 
336  template < typename LhsT, typename RhsT >
337  Derived& operator= ( const Addition< LhsT, RhsT > &add ) ;
338 
339  template < typename Expression >
340  Derived& operator= ( const NarySum< Expression > &sum ) ;
341 
342  template < typename OtherDerived >
343  Derived& operator= ( const Scaling< OtherDerived > &scaling )
344  {
345  Evaluator< typename Scaling< OtherDerived >::Operand::ObjectType > rhs ( scaling.operand.object ) ;
346  return assign< Scaling< OtherDerived >::transposeOperand >( *rhs, scaling.operand.scaling ) ;
347  }
348 
350  template< typename OtherDerived >
351  void cloneDimensions( const BlockMatrixBase< OtherDerived > &source ) ;
352 
354  template< typename OtherDerived >
355  void cloneStructure( const SparseBlockMatrixBase< OtherDerived > &source ) ;
356 
359  void cloneIndex( const MajorIndexType &index ) ;
360 
362 
364 
365 
366  ConstTransposeReturnType transpose() const { return Transpose< Derived >( derived() ) ; }
367 
368  template < bool Transpose, typename RhsT, typename ResT >
369  void multiply( const RhsT& rhs, ResT& res, Scalar alpha = 1, Scalar beta = 0 ) const ;
370 
371  template < typename RhsT, typename ResT >
372  void splitRowMultiply( const Index row, const RhsT& rhs, ResT& res ) const ;
373 
374  template < bool DoTranspose, typename RhsT, typename ResT, typename PreOp >
375  void rowMultiplyPrecompose( const Index row, const RhsT& rhs, ResT& res, const PreOp &op ) const ;
376 
377  template < bool DoTranspose, typename RhsT, typename ResT, typename PostOp >
378  void colMultiplyPostcompose( const Index row, const RhsT& rhs, ResT& res, const PostOp &op ) const ;
379 
380  template < bool ColWise, typename LhsT, typename RhsT >
381  void setFromProduct( const Product< LhsT, RhsT > &prod ) ;
382 
384  Derived& scale( Scalar alpha ) ;
385 
387  template < bool Transpose, typename OtherDerived >
388  Derived& add( const SparseBlockMatrixBase< OtherDerived > &rhs, Scalar alpha = 1) ;
389 
391  Derived& operator *= ( Scalar alpha ) { return scale( alpha ) ; }
393  Derived& operator /= ( Scalar alpha ) { return scale( 1./alpha ) ; }
394 
397  { return Scaling< Derived >( derived(), -1 ) ; }
398 
400  template < typename OtherDerived >
402  {
403  Evaluator< OtherDerived > rhs ( source.derived() ) ;
404  return add< OtherDerived::is_transposed >( *rhs ) ;
405  }
406 
408  template < typename OtherDerived >
410  {
411  Evaluator< OtherDerived > rhs ( source.derived() ) ;
412  return add< OtherDerived::is_transposed >( *rhs, -1 ) ;
413  }
414 
416 
418 
419 
420 #ifdef BOGUS_WITH_BOOST_SERIALIZATION
421  template < typename Archive >
422  void serialize( Archive & ar, const unsigned int file_version ) ;
423 #endif
424 
426 
428 
429 
431  Derived& prune( const Scalar precision = 0 ) ;
432 
434 
440  Derived& applyPermutation( const std::size_t* indices ) ;
441 
443 
445 
446  CopyResultType Zero() const
447  {
448  CopyResultType m ;
449  m.cloneDimensions( *this ) ;
450  return m.setZero() ;
451  }
452 
453  CopyResultType Identity() const
454  {
455  CopyResultType m ;
456  m.cloneDimensions( *this ) ;
457  return m.setIdentity() ;
458  }
460 
462 
463 
465 
467  MajorIndexType& majorIndex() { return m_majorIndex ; }
468 
470  void prealloc( std::size_t nBlocks ) ;
471 
473  const TransposeArrayType& transposeBlocks() const
474  { return m_transposeBlocks ; }
475 
477  const TransposeBlockType* transposeData() const
478  { return data_pointer(m_transposeBlocks) ; }
479 
481  const Index* rowOffsets() const { return colMajorIndex().innerOffsetsData() ; }
482 
484  const Index* colOffsets() const { return rowMajorIndex().innerOffsetsData() ; }
485 
487  const Lock& lock() const { return m_lock ; }
488 
490 
491 protected:
492 
493  enum {
494  is_bsr_compatible =
495  Base::has_fixed_size_blocks && Base::has_square_or_dynamic_blocks &&
496  Traits::is_compressed && !Traits::is_col_major &&
498  } ;
499 
500  using Base::m_cols ;
501  using Base::m_rows ;
502  using Base::m_blocks ;
503 
504  typedef SparseBlockMatrixFinalizer< Traits::is_symmetric > Finalizer ;
505  friend struct SparseBlockIndexGetter< Derived, true > ;
506  friend struct SparseBlockIndexGetter< Derived, false > ;
507 
508  SparseBlockMatrixBase() ;
509 
511  template< bool EnforceThreadSafety >
512  BlockType& allocateBlock( BlockPtr &ptr ) ;
513 
514  void computeMinorIndex( MinorIndexType &cmIndex) const ;
515 
516  const MinorIndexType& getOrComputeMinorIndex( MinorIndexType &tempIndex) const ;
517 
518  ColIndexType& colMajorIndex() ;
519  RowIndexType& rowMajorIndex() ;
520 
521  template< typename IndexT >
522  void setInnerOffets( IndexT& index, const Index nBlocks, const unsigned* blockSizes ) const ;
523 
524  TransposeBlockType* transposeData()
525  { return data_pointer(m_transposeBlocks) ; }
526 
527 
528  TransposeArrayType m_transposeBlocks ;
529 
530  MajorIndexType m_majorIndex ;
531 
532  MinorIndexType m_minorIndex ;
533 
534  TransposeIndexType m_transposeIndex ;
535 
536  Lock m_lock ;
537 } ;
538 
539 }
540 
541 #endif // SPARSEBLOCKMATRIX_HH
const Traits::BlocksArrayType & blocks() const
Access to blocks data.
Definition: BlockMatrixBase.hpp:89
Base class for dense and sparse block matrices, thought dense don&#39;t exist yet.
Definition: BlockMatrixBase.hpp:22
BlockType & allocateBlock(BlockPtr &ptr)
Pushes a block at the back of m_blocks.
const Lock & lock() const
Reference to matrix private mutex.
Definition: SparseBlockMatrixBase.hpp:487
std::size_t nBlocks() const
Returns the number of (non-zero) blocks of the matrix.
Definition: SparseBlockMatrixBase.hpp:240
void cloneIndex(const MajorIndexType &index)
Definition: Traits.hpp:19
Uncompressed sparse block index.
Definition: SparseBlockIndex.hpp:86
Derived & operator*=(Scalar alpha)
Coeff-wise multiplication with a scalar.
Definition: SparseBlockMatrixBase.hpp:391
void prealloc(std::size_t nBlocks)
Resizes m_blocks.
Derived & scale(Scalar alpha)
Performs *this *= alpha.
const TransposeBlockType * transposeData() const
Returns the blocks that have been created by cacheTranspose(), as a raw pointer.
Definition: SparseBlockMatrixBase.hpp:477
void setCols(const Index nBlocks, const Index colsPerBlock)
Same, setting each block to have exactly rowsPerBlock.
Definition: SparseBlockMatrixBase.hpp:148
Definition: SparseBlockMatrixBase.hpp:29
const Derived & derived() const
Returns a const reference to the implementation.
Definition: BlockObjectBase.hpp:25
BlockType & insertBack(Index row, Index col)
Inserts a block at the end of the matrix, and returns a reference to it.
Definition: SparseBlockMatrixBase.hpp:183
BlockType & insertAndResize(Index row, Index col)
Convenience method that insert() a block and immediately resize it according to the dimensions given ...
void reserve(std::size_t nBlocks)
Reserve enouch memory to accomodate nBlocks.
Definition: SparseBlockMatrixBase.hpp:171
void setCols(const Index nBlocks)
Same, deducing the (constant) number of rows per block from the BlockType.
Definition: SparseBlockMatrixBase.hpp:153
Derived & operator+=(const BlockObjectBase< OtherDerived > &source)
Adds another SparseBlockMatrixBase to this one.
Definition: SparseBlockMatrixBase.hpp:401
Definition: SparseBlockIndexComputer.hpp:23
Sparse Block Matrix.
Definition: SparseBlockMatrix.hpp:57
Index cols() const
Returns the total number of columns of the matrix ( expanding blocks )
Definition: BlockMatrixBase.hpp:83
Base class for matrix-like objects that define a block structure, but not a block type...
Definition: IterableBlockObject.hpp:22
void clear()
Clears the matrix.
const BlockType & block(BlockPtr ptr) const
Returns a reference to a block using the result from blockPtr() or diagonalBlockPtr() ...
Definition: BlockMatrixBase.hpp:70
BlockType & insertBackAndResize(Index row, Index col)
Convenience method that insertBack() a block and immediately resize it according to the dimensions gi...
bool computeMinorIndex()
Computes the minor index of the matrix.
Definition: Expressions.hpp:268
Representation of the null matrix.
Definition: Zero.hpp:21
BlockPtr blockPtr(Index row, Index col) const
Return a BlockPtr to the block a (row, col) or InvalidBlockPtr if it does not exist.
BlockType & diagonal(const Index row)
Definition: BlockMatrixBase.hpp:96
Definition: SparseBlockMatrixBase.hpp:70
static const BlockPtr InvalidBlockPtr
Return value of blockPtr( Index, Index ) for non-existing block.
Definition: BlockMatrixBase.hpp:38
Index rows() const
Returns the total number of rows of the matrix ( expanding blocks )
Definition: BlockMatrixBase.hpp:81
bool transposeCached() const
Returns whether the transpose has been cached.
Definition: SparseBlockMatrixBase.hpp:283
Definition: Traits.hpp:29
bool empty() const
Returns whether the matrix is empty.
Definition: SparseBlockMatrixBase.hpp:244
Derived & operator-=(const BlockObjectBase< OtherDerived > &source)
Substracts another SparseBlockMatrixBase from this one.
Definition: SparseBlockMatrixBase.hpp:409
bool valid
Whether this index is currently valid.
Definition: SparseBlockIndex.hpp:40
Evaluates an expression inside a temporary if necessary, otherwise returns directly a matrix referenc...
Definition: Evaluators.hpp:20
Definition: Lock.hpp:51
MajorIndexType & majorIndex()
Direct access to major index.
Definition: SparseBlockMatrixBase.hpp:467
void cloneDimensions(const BlockMatrixBase< OtherDerived > &source)
Clones the dimensions ( number of rows/cols blocks and rows/cols per block ) of source.
void setRows(const Index nBlocks, const Index rowsPerBlock)
Same, setting each block to have exactly rowsPerBlock.
Definition: SparseBlockMatrixBase.hpp:126
Derived & setZero()
Definition: SparseBlockMatrixBase.hpp:233
Derived & applyPermutation(const std::size_t *indices)
Sets *this = P * (*this) * P.transpose().
Derived & setIdentity()
Calls set_identity() on each diagonal block, discard off-diagonal blocks.
void setCols(const Index nBlocks, const unsigned *colsPerBlock)
Defines the column structure of the matrix.
Base class for anything block.
Definition: BlockObjectBase.hpp:22
void eachBlockOf(const Index outer, Func func) const
Iterates over each block of a given row or col. Calls func( inner, block )
Default container type, that should resizable and use contiguous storage.
Definition: Traits.hpp:23
Derived & add(const SparseBlockMatrixBase< OtherDerived > &rhs, Scalar alpha=1)
Performs *this += alpha * rhs (SAXPY)
Derived & setBlocksToZero()
Sets all allocated blocks to zero. Does not update index.
const TransposeArrayType & transposeBlocks() const
Returns the blocks that have been created by cacheTranspose()
Definition: SparseBlockMatrixBase.hpp:473
Store and index blocks in a column major way.
Definition: Constants.hpp:58
void setRows(const Index nBlocks, const unsigned *rowsPerBlock)
Defines the row structure of the matrix.
Derived & operator/=(Scalar alpha)
Coeff-wise division with a scalar.
Definition: SparseBlockMatrixBase.hpp:393
void setRows(const Index nBlocks)
Same, deducing the (constant) number of rows per block from the BlockType.
Definition: SparseBlockMatrixBase.hpp:131
BlockType & insertByOuterInner(Index outer, Index inner)
Insert a block, specifying directily the outer and inner indices instead of row and column...
BlockType & insert(Index row, Index col)
Insert a block anywhere in the matrix, and returns a reference to it.
Definition: SparseBlockMatrixBase.hpp:208
const Index * rowOffsets() const
Returns an array containing the first index of each row.
Definition: SparseBlockMatrixBase.hpp:481
BlockPtr diagonalBlockPtr(Index row) const
Return a BlockPtr to the block a (row, row) or InvalidBlockPtr if it does not exist.
void setRows(const std::vector< unsigned > &rowsPerBlock)
Same, using a std::vector.
Definition: SparseBlockMatrixBase.hpp:121
void cacheTranspose()
Computes and caches the tranpose of the matrix.
const Index * colOffsets() const
Returns an array containing the first index of each column.
Definition: SparseBlockMatrixBase.hpp:484
BlockType TransposeStorageType
Type for storing the result of transpose_block( BlockType ), useful for cacheTranspose() ...
Definition: Traits.hpp:33
Derived & assign(const SparseBlockMatrixBase< OtherDerived > &source, const Scalar scale=1)
Performs ( *this = scale * source ) or ( *this = scale * source.transpose() )
void cloneStructure(const SparseBlockMatrixBase< OtherDerived > &source)
Clones the dimensions and the indexes of source.
Scaling< Derived > operator-() const
Unary minus.
Definition: SparseBlockMatrixBase.hpp:396
Base class for SparseBlockMatrix.
Definition: SparseBlockMatrixBase.hpp:36
InnerIterator innerIterator(Index outer) const
Returns an InnerIterator to efficiently browse matrix.
Definition: SparseBlockMatrixBase.hpp:295
void setCols(const std::vector< unsigned > &colsPerBlock)
Same, using a std::vector.
Definition: SparseBlockMatrixBase.hpp:143
Derived & prune(const Scalar precision=0)
Removes all blocks for which is_zero( block, precision ) is true.