1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud (at) inria.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_SPARSE_BLOCK_H 11 #define EIGEN_SPARSE_BLOCK_H 12 13 namespace Eigen { 14 15 template<typename XprType, int BlockRows, int BlockCols> 16 class BlockImpl<XprType,BlockRows,BlockCols,true,Sparse> 17 : public SparseMatrixBase<Block<XprType,BlockRows,BlockCols,true> > 18 { 19 typedef typename internal::remove_all<typename XprType::Nested>::type _MatrixTypeNested; 20 typedef Block<XprType, BlockRows, BlockCols, true> BlockType; 21 public: 22 enum { IsRowMajor = internal::traits<BlockType>::IsRowMajor }; 23 protected: 24 enum { OuterSize = IsRowMajor ? BlockRows : BlockCols }; 25 public: 26 EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType) 27 28 class InnerIterator: public XprType::InnerIterator 29 { 30 typedef typename BlockImpl::Index Index; 31 public: 32 inline InnerIterator(const BlockType& xpr, Index outer) 33 : XprType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) 34 {} 35 inline Index row() const { return IsRowMajor ? m_outer : this->index(); } 36 inline Index col() const { return IsRowMajor ? this->index() : m_outer; } 37 protected: 38 Index m_outer; 39 }; 40 class ReverseInnerIterator: public XprType::ReverseInnerIterator 41 { 42 typedef typename BlockImpl::Index Index; 43 public: 44 inline ReverseInnerIterator(const BlockType& xpr, Index outer) 45 : XprType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) 46 {} 47 inline Index row() const { return IsRowMajor ? m_outer : this->index(); } 48 inline Index col() const { return IsRowMajor ? this->index() : m_outer; } 49 protected: 50 Index m_outer; 51 }; 52 53 inline BlockImpl(const XprType& xpr, int i) 54 : m_matrix(xpr), m_outerStart(i), m_outerSize(OuterSize) 55 {} 56 57 inline BlockImpl(const XprType& xpr, int startRow, int startCol, int blockRows, int blockCols) 58 : m_matrix(xpr), m_outerStart(IsRowMajor ? startRow : startCol), m_outerSize(IsRowMajor ? blockRows : blockCols) 59 {} 60 61 EIGEN_STRONG_INLINE Index rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); } 62 EIGEN_STRONG_INLINE Index cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); } 63 64 protected: 65 66 typename XprType::Nested m_matrix; 67 Index m_outerStart; 68 const internal::variable_if_dynamic<Index, OuterSize> m_outerSize; 69 70 EIGEN_INHERIT_ASSIGNMENT_OPERATORS(BlockImpl) 71 }; 72 73 74 /*************************************************************************** 75 * specialisation for SparseMatrix 76 ***************************************************************************/ 77 78 template<typename _Scalar, int _Options, typename _Index, int BlockRows, int BlockCols> 79 class BlockImpl<SparseMatrix<_Scalar, _Options, _Index>,BlockRows,BlockCols,true,Sparse> 80 : public SparseMatrixBase<Block<SparseMatrix<_Scalar, _Options, _Index>,BlockRows,BlockCols,true> > 81 { 82 typedef SparseMatrix<_Scalar, _Options, _Index> SparseMatrixType; 83 typedef typename internal::remove_all<typename SparseMatrixType::Nested>::type _MatrixTypeNested; 84 typedef Block<SparseMatrixType, BlockRows, BlockCols, true> BlockType; 85 public: 86 enum { IsRowMajor = internal::traits<BlockType>::IsRowMajor }; 87 EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType) 88 protected: 89 enum { OuterSize = IsRowMajor ? BlockRows : BlockCols }; 90 public: 91 92 class InnerIterator: public SparseMatrixType::InnerIterator 93 { 94 public: 95 inline InnerIterator(const BlockType& xpr, Index outer) 96 : SparseMatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) 97 {} 98 inline Index row() const { return IsRowMajor ? m_outer : this->index(); } 99 inline Index col() const { return IsRowMajor ? this->index() : m_outer; } 100 protected: 101 Index m_outer; 102 }; 103 class ReverseInnerIterator: public SparseMatrixType::ReverseInnerIterator 104 { 105 public: 106 inline ReverseInnerIterator(const BlockType& xpr, Index outer) 107 : SparseMatrixType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) 108 {} 109 inline Index row() const { return IsRowMajor ? m_outer : this->index(); } 110 inline Index col() const { return IsRowMajor ? this->index() : m_outer; } 111 protected: 112 Index m_outer; 113 }; 114 115 inline BlockImpl(const SparseMatrixType& xpr, int i) 116 : m_matrix(xpr), m_outerStart(i), m_outerSize(OuterSize) 117 {} 118 119 inline BlockImpl(const SparseMatrixType& xpr, int startRow, int startCol, int blockRows, int blockCols) 120 : m_matrix(xpr), m_outerStart(IsRowMajor ? startRow : startCol), m_outerSize(IsRowMajor ? blockRows : blockCols) 121 {} 122 123 template<typename OtherDerived> 124 inline BlockType& operator=(const SparseMatrixBase<OtherDerived>& other) 125 { 126 typedef typename internal::remove_all<typename SparseMatrixType::Nested>::type _NestedMatrixType; 127 _NestedMatrixType& matrix = const_cast<_NestedMatrixType&>(m_matrix);; 128 // This assignement is slow if this vector set is not empty 129 // and/or it is not at the end of the nonzeros of the underlying matrix. 130 131 // 1 - eval to a temporary to avoid transposition and/or aliasing issues 132 SparseMatrix<Scalar, IsRowMajor ? RowMajor : ColMajor, Index> tmp(other); 133 134 // 2 - let's check whether there is enough allocated memory 135 Index nnz = tmp.nonZeros(); 136 Index start = m_outerStart==0 ? 0 : matrix.outerIndexPtr()[m_outerStart]; // starting position of the current block 137 Index end = m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()]; // ending posiiton of the current block 138 Index block_size = end - start; // available room in the current block 139 Index tail_size = m_matrix.outerIndexPtr()[m_matrix.outerSize()] - end; 140 141 Index free_size = m_matrix.isCompressed() 142 ? Index(matrix.data().allocatedSize()) + block_size 143 : block_size; 144 145 if(nnz>free_size) 146 { 147 // realloc manually to reduce copies 148 typename SparseMatrixType::Storage newdata(m_matrix.data().allocatedSize() - block_size + nnz); 149 150 std::memcpy(&newdata.value(0), &m_matrix.data().value(0), start*sizeof(Scalar)); 151 std::memcpy(&newdata.index(0), &m_matrix.data().index(0), start*sizeof(Index)); 152 153 std::memcpy(&newdata.value(start), &tmp.data().value(0), nnz*sizeof(Scalar)); 154 std::memcpy(&newdata.index(start), &tmp.data().index(0), nnz*sizeof(Index)); 155 156 std::memcpy(&newdata.value(start+nnz), &matrix.data().value(end), tail_size*sizeof(Scalar)); 157 std::memcpy(&newdata.index(start+nnz), &matrix.data().index(end), tail_size*sizeof(Index)); 158 159 newdata.resize(m_matrix.outerIndexPtr()[m_matrix.outerSize()] - block_size + nnz); 160 161 matrix.data().swap(newdata); 162 } 163 else 164 { 165 // no need to realloc, simply copy the tail at its respective position and insert tmp 166 matrix.data().resize(start + nnz + tail_size); 167 168 std::memmove(&matrix.data().value(start+nnz), &matrix.data().value(end), tail_size*sizeof(Scalar)); 169 std::memmove(&matrix.data().index(start+nnz), &matrix.data().index(end), tail_size*sizeof(Index)); 170 171 std::memcpy(&matrix.data().value(start), &tmp.data().value(0), nnz*sizeof(Scalar)); 172 std::memcpy(&matrix.data().index(start), &tmp.data().index(0), nnz*sizeof(Index)); 173 } 174 175 // update innerNonZeros 176 if(!m_matrix.isCompressed()) 177 for(Index j=0; j<m_outerSize.value(); ++j) 178 matrix.innerNonZeroPtr()[m_outerStart+j] = tmp.innerVector(j).nonZeros(); 179 180 // update outer index pointers 181 Index p = start; 182 for(Index k=0; k<m_outerSize.value(); ++k) 183 { 184 matrix.outerIndexPtr()[m_outerStart+k] = p; 185 p += tmp.innerVector(k).nonZeros(); 186 } 187 std::ptrdiff_t offset = nnz - block_size; 188 for(Index k = m_outerStart + m_outerSize.value(); k<=matrix.outerSize(); ++k) 189 { 190 matrix.outerIndexPtr()[k] += offset; 191 } 192 193 return derived(); 194 } 195 196 inline BlockType& operator=(const BlockType& other) 197 { 198 return operator=<BlockType>(other); 199 } 200 201 inline const Scalar* valuePtr() const 202 { return m_matrix.valuePtr() + m_matrix.outerIndexPtr()[m_outerStart]; } 203 inline Scalar* valuePtr() 204 { return m_matrix.const_cast_derived().valuePtr() + m_matrix.outerIndexPtr()[m_outerStart]; } 205 206 inline const Index* innerIndexPtr() const 207 { return m_matrix.innerIndexPtr() + m_matrix.outerIndexPtr()[m_outerStart]; } 208 inline Index* innerIndexPtr() 209 { return m_matrix.const_cast_derived().innerIndexPtr() + m_matrix.outerIndexPtr()[m_outerStart]; } 210 211 inline const Index* outerIndexPtr() const 212 { return m_matrix.outerIndexPtr() + m_outerStart; } 213 inline Index* outerIndexPtr() 214 { return m_matrix.const_cast_derived().outerIndexPtr() + m_outerStart; } 215 216 Index nonZeros() const 217 { 218 if(m_matrix.isCompressed()) 219 return std::size_t(m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()]) 220 - std::size_t(m_matrix.outerIndexPtr()[m_outerStart]); 221 else if(m_outerSize.value()==0) 222 return 0; 223 else 224 return Map<const Matrix<Index,OuterSize,1> >(m_matrix.innerNonZeroPtr()+m_outerStart, m_outerSize.value()).sum(); 225 } 226 227 const Scalar& lastCoeff() const 228 { 229 EIGEN_STATIC_ASSERT_VECTOR_ONLY(BlockImpl); 230 eigen_assert(nonZeros()>0); 231 if(m_matrix.isCompressed()) 232 return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart+1]-1]; 233 else 234 return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart]+m_matrix.innerNonZeroPtr()[m_outerStart]-1]; 235 } 236 237 EIGEN_STRONG_INLINE Index rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); } 238 EIGEN_STRONG_INLINE Index cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); } 239 240 protected: 241 242 typename SparseMatrixType::Nested m_matrix; 243 Index m_outerStart; 244 const internal::variable_if_dynamic<Index, OuterSize> m_outerSize; 245 246 }; 247 248 //---------- 249 250 /** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this 251 * is col-major (resp. row-major). 252 */ 253 template<typename Derived> 254 typename SparseMatrixBase<Derived>::InnerVectorReturnType SparseMatrixBase<Derived>::innerVector(Index outer) 255 { return InnerVectorReturnType(derived(), outer); } 256 257 /** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this 258 * is col-major (resp. row-major). Read-only. 259 */ 260 template<typename Derived> 261 const typename SparseMatrixBase<Derived>::ConstInnerVectorReturnType SparseMatrixBase<Derived>::innerVector(Index outer) const 262 { return ConstInnerVectorReturnType(derived(), outer); } 263 264 /** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this 265 * is col-major (resp. row-major). 266 */ 267 template<typename Derived> 268 Block<Derived,Dynamic,Dynamic,true> SparseMatrixBase<Derived>::innerVectors(Index outerStart, Index outerSize) 269 { 270 return Block<Derived,Dynamic,Dynamic,true>(derived(), 271 IsRowMajor ? outerStart : 0, IsRowMajor ? 0 : outerStart, 272 IsRowMajor ? outerSize : rows(), IsRowMajor ? cols() : outerSize); 273 274 } 275 276 /** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this 277 * is col-major (resp. row-major). Read-only. 278 */ 279 template<typename Derived> 280 const Block<const Derived,Dynamic,Dynamic,true> SparseMatrixBase<Derived>::innerVectors(Index outerStart, Index outerSize) const 281 { 282 return Block<const Derived,Dynamic,Dynamic,true>(derived(), 283 IsRowMajor ? outerStart : 0, IsRowMajor ? 0 : outerStart, 284 IsRowMajor ? outerSize : rows(), IsRowMajor ? cols() : outerSize); 285 286 } 287 288 /** Generic implementation of sparse Block expression. 289 * Real-only. 290 */ 291 template<typename XprType, int BlockRows, int BlockCols, bool InnerPanel> 292 class BlockImpl<XprType,BlockRows,BlockCols,InnerPanel,Sparse> 293 : public SparseMatrixBase<Block<XprType,BlockRows,BlockCols,InnerPanel> >, internal::no_assignment_operator 294 { 295 typedef typename internal::remove_all<typename XprType::Nested>::type _MatrixTypeNested; 296 typedef Block<XprType, BlockRows, BlockCols, InnerPanel> BlockType; 297 public: 298 enum { IsRowMajor = internal::traits<BlockType>::IsRowMajor }; 299 EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType) 300 301 /** Column or Row constructor 302 */ 303 inline BlockImpl(const XprType& xpr, int i) 304 : m_matrix(xpr), 305 m_startRow( (BlockRows==1) && (BlockCols==XprType::ColsAtCompileTime) ? i : 0), 306 m_startCol( (BlockRows==XprType::RowsAtCompileTime) && (BlockCols==1) ? i : 0), 307 m_blockRows(xpr.rows()), 308 m_blockCols(xpr.cols()) 309 {} 310 311 /** Dynamic-size constructor 312 */ 313 inline BlockImpl(const XprType& xpr, int startRow, int startCol, int blockRows, int blockCols) 314 : m_matrix(xpr), m_startRow(startRow), m_startCol(startCol), m_blockRows(blockRows), m_blockCols(blockCols) 315 {} 316 317 inline int rows() const { return m_blockRows.value(); } 318 inline int cols() const { return m_blockCols.value(); } 319 320 inline Scalar& coeffRef(int row, int col) 321 { 322 return m_matrix.const_cast_derived() 323 .coeffRef(row + m_startRow.value(), col + m_startCol.value()); 324 } 325 326 inline const Scalar coeff(int row, int col) const 327 { 328 return m_matrix.coeff(row + m_startRow.value(), col + m_startCol.value()); 329 } 330 331 inline Scalar& coeffRef(int index) 332 { 333 return m_matrix.const_cast_derived() 334 .coeffRef(m_startRow.value() + (RowsAtCompileTime == 1 ? 0 : index), 335 m_startCol.value() + (RowsAtCompileTime == 1 ? index : 0)); 336 } 337 338 inline const Scalar coeff(int index) const 339 { 340 return m_matrix 341 .coeff(m_startRow.value() + (RowsAtCompileTime == 1 ? 0 : index), 342 m_startCol.value() + (RowsAtCompileTime == 1 ? index : 0)); 343 } 344 345 inline const _MatrixTypeNested& nestedExpression() const { return m_matrix; } 346 347 class InnerIterator : public _MatrixTypeNested::InnerIterator 348 { 349 typedef typename _MatrixTypeNested::InnerIterator Base; 350 const BlockType& m_block; 351 Index m_end; 352 public: 353 354 EIGEN_STRONG_INLINE InnerIterator(const BlockType& block, Index outer) 355 : Base(block.derived().nestedExpression(), outer + (IsRowMajor ? block.m_startRow.value() : block.m_startCol.value())), 356 m_block(block), 357 m_end(IsRowMajor ? block.m_startCol.value()+block.m_blockCols.value() : block.m_startRow.value()+block.m_blockRows.value()) 358 { 359 while( (Base::operator bool()) && (Base::index() < (IsRowMajor ? m_block.m_startCol.value() : m_block.m_startRow.value())) ) 360 Base::operator++(); 361 } 362 363 inline Index index() const { return Base::index() - (IsRowMajor ? m_block.m_startCol.value() : m_block.m_startRow.value()); } 364 inline Index outer() const { return Base::outer() - (IsRowMajor ? m_block.m_startRow.value() : m_block.m_startCol.value()); } 365 inline Index row() const { return Base::row() - m_block.m_startRow.value(); } 366 inline Index col() const { return Base::col() - m_block.m_startCol.value(); } 367 368 inline operator bool() const { return Base::operator bool() && Base::index() < m_end; } 369 }; 370 class ReverseInnerIterator : public _MatrixTypeNested::ReverseInnerIterator 371 { 372 typedef typename _MatrixTypeNested::ReverseInnerIterator Base; 373 const BlockType& m_block; 374 Index m_begin; 375 public: 376 377 EIGEN_STRONG_INLINE ReverseInnerIterator(const BlockType& block, Index outer) 378 : Base(block.derived().nestedExpression(), outer + (IsRowMajor ? block.m_startRow.value() : block.m_startCol.value())), 379 m_block(block), 380 m_begin(IsRowMajor ? block.m_startCol.value() : block.m_startRow.value()) 381 { 382 while( (Base::operator bool()) && (Base::index() >= (IsRowMajor ? m_block.m_startCol.value()+block.m_blockCols.value() : m_block.m_startRow.value()+block.m_blockRows.value())) ) 383 Base::operator--(); 384 } 385 386 inline Index index() const { return Base::index() - (IsRowMajor ? m_block.m_startCol.value() : m_block.m_startRow.value()); } 387 inline Index outer() const { return Base::outer() - (IsRowMajor ? m_block.m_startRow.value() : m_block.m_startCol.value()); } 388 inline Index row() const { return Base::row() - m_block.m_startRow.value(); } 389 inline Index col() const { return Base::col() - m_block.m_startCol.value(); } 390 391 inline operator bool() const { return Base::operator bool() && Base::index() >= m_begin; } 392 }; 393 protected: 394 friend class InnerIterator; 395 friend class ReverseInnerIterator; 396 397 EIGEN_INHERIT_ASSIGNMENT_OPERATORS(BlockImpl) 398 399 typename XprType::Nested m_matrix; 400 const internal::variable_if_dynamic<Index, XprType::RowsAtCompileTime == 1 ? 0 : Dynamic> m_startRow; 401 const internal::variable_if_dynamic<Index, XprType::ColsAtCompileTime == 1 ? 0 : Dynamic> m_startCol; 402 const internal::variable_if_dynamic<Index, RowsAtCompileTime> m_blockRows; 403 const internal::variable_if_dynamic<Index, ColsAtCompileTime> m_blockCols; 404 405 }; 406 407 } // end namespace Eigen 408 409 #endif // EIGEN_SPARSE_BLOCK_H 410