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@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 13namespace Eigen { 14 15template<typename XprType, int BlockRows, int BlockCols> 16class 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; 21public: 22 enum { IsRowMajor = internal::traits<BlockType>::IsRowMajor }; 23protected: 24 enum { OuterSize = IsRowMajor ? BlockRows : BlockCols }; 25public: 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 private: 72 Index nonZeros() const; 73}; 74 75 76/*************************************************************************** 77* specialisation for SparseMatrix 78***************************************************************************/ 79 80template<typename _Scalar, int _Options, typename _Index, int BlockRows, int BlockCols> 81class BlockImpl<SparseMatrix<_Scalar, _Options, _Index>,BlockRows,BlockCols,true,Sparse> 82 : public SparseMatrixBase<Block<SparseMatrix<_Scalar, _Options, _Index>,BlockRows,BlockCols,true> > 83{ 84 typedef SparseMatrix<_Scalar, _Options, _Index> SparseMatrixType; 85 typedef typename internal::remove_all<typename SparseMatrixType::Nested>::type _MatrixTypeNested; 86 typedef Block<SparseMatrixType, BlockRows, BlockCols, true> BlockType; 87 typedef Block<const SparseMatrixType, BlockRows, BlockCols, true> ConstBlockType; 88public: 89 enum { IsRowMajor = internal::traits<BlockType>::IsRowMajor }; 90 EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType) 91protected: 92 enum { OuterSize = IsRowMajor ? BlockRows : BlockCols }; 93public: 94 95 class InnerIterator: public SparseMatrixType::InnerIterator 96 { 97 public: 98 inline InnerIterator(const BlockType& xpr, Index outer) 99 : SparseMatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) 100 {} 101 inline Index row() const { return IsRowMajor ? m_outer : this->index(); } 102 inline Index col() const { return IsRowMajor ? this->index() : m_outer; } 103 protected: 104 Index m_outer; 105 }; 106 class ReverseInnerIterator: public SparseMatrixType::ReverseInnerIterator 107 { 108 public: 109 inline ReverseInnerIterator(const BlockType& xpr, Index outer) 110 : SparseMatrixType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) 111 {} 112 inline Index row() const { return IsRowMajor ? m_outer : this->index(); } 113 inline Index col() const { return IsRowMajor ? this->index() : m_outer; } 114 protected: 115 Index m_outer; 116 }; 117 118 inline BlockImpl(const SparseMatrixType& xpr, int i) 119 : m_matrix(xpr), m_outerStart(i), m_outerSize(OuterSize) 120 {} 121 122 inline BlockImpl(const SparseMatrixType& xpr, int startRow, int startCol, int blockRows, int blockCols) 123 : m_matrix(xpr), m_outerStart(IsRowMajor ? startRow : startCol), m_outerSize(IsRowMajor ? blockRows : blockCols) 124 {} 125 126 template<typename OtherDerived> 127 inline BlockType& operator=(const SparseMatrixBase<OtherDerived>& other) 128 { 129 typedef typename internal::remove_all<typename SparseMatrixType::Nested>::type _NestedMatrixType; 130 _NestedMatrixType& matrix = const_cast<_NestedMatrixType&>(m_matrix);; 131 // This assignement is slow if this vector set is not empty 132 // and/or it is not at the end of the nonzeros of the underlying matrix. 133 134 // 1 - eval to a temporary to avoid transposition and/or aliasing issues 135 SparseMatrix<Scalar, IsRowMajor ? RowMajor : ColMajor, Index> tmp(other); 136 137 // 2 - let's check whether there is enough allocated memory 138 Index nnz = tmp.nonZeros(); 139 Index start = m_outerStart==0 ? 0 : matrix.outerIndexPtr()[m_outerStart]; // starting position of the current block 140 Index end = m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()]; // ending posiiton of the current block 141 Index block_size = end - start; // available room in the current block 142 Index tail_size = m_matrix.outerIndexPtr()[m_matrix.outerSize()] - end; 143 144 Index free_size = m_matrix.isCompressed() 145 ? Index(matrix.data().allocatedSize()) + block_size 146 : block_size; 147 148 if(nnz>free_size) 149 { 150 // realloc manually to reduce copies 151 typename SparseMatrixType::Storage newdata(m_matrix.data().allocatedSize() - block_size + nnz); 152 153 std::memcpy(&newdata.value(0), &m_matrix.data().value(0), start*sizeof(Scalar)); 154 std::memcpy(&newdata.index(0), &m_matrix.data().index(0), start*sizeof(Index)); 155 156 std::memcpy(&newdata.value(start), &tmp.data().value(0), nnz*sizeof(Scalar)); 157 std::memcpy(&newdata.index(start), &tmp.data().index(0), nnz*sizeof(Index)); 158 159 std::memcpy(&newdata.value(start+nnz), &matrix.data().value(end), tail_size*sizeof(Scalar)); 160 std::memcpy(&newdata.index(start+nnz), &matrix.data().index(end), tail_size*sizeof(Index)); 161 162 newdata.resize(m_matrix.outerIndexPtr()[m_matrix.outerSize()] - block_size + nnz); 163 164 matrix.data().swap(newdata); 165 } 166 else 167 { 168 // no need to realloc, simply copy the tail at its respective position and insert tmp 169 matrix.data().resize(start + nnz + tail_size); 170 171 std::memmove(&matrix.data().value(start+nnz), &matrix.data().value(end), tail_size*sizeof(Scalar)); 172 std::memmove(&matrix.data().index(start+nnz), &matrix.data().index(end), tail_size*sizeof(Index)); 173 174 std::memcpy(&matrix.data().value(start), &tmp.data().value(0), nnz*sizeof(Scalar)); 175 std::memcpy(&matrix.data().index(start), &tmp.data().index(0), nnz*sizeof(Index)); 176 } 177 178 // update innerNonZeros 179 if(!m_matrix.isCompressed()) 180 for(Index j=0; j<m_outerSize.value(); ++j) 181 matrix.innerNonZeroPtr()[m_outerStart+j] = tmp.innerVector(j).nonZeros(); 182 183 // update outer index pointers 184 Index p = start; 185 for(Index k=0; k<m_outerSize.value(); ++k) 186 { 187 matrix.outerIndexPtr()[m_outerStart+k] = p; 188 p += tmp.innerVector(k).nonZeros(); 189 } 190 std::ptrdiff_t offset = nnz - block_size; 191 for(Index k = m_outerStart + m_outerSize.value(); k<=matrix.outerSize(); ++k) 192 { 193 matrix.outerIndexPtr()[k] += offset; 194 } 195 196 return derived(); 197 } 198 199 inline BlockType& operator=(const BlockType& other) 200 { 201 return operator=<BlockType>(other); 202 } 203 204 inline const Scalar* valuePtr() const 205 { return m_matrix.valuePtr() + m_matrix.outerIndexPtr()[m_outerStart]; } 206 inline Scalar* valuePtr() 207 { return m_matrix.const_cast_derived().valuePtr() + m_matrix.outerIndexPtr()[m_outerStart]; } 208 209 inline const Index* innerIndexPtr() const 210 { return m_matrix.innerIndexPtr() + m_matrix.outerIndexPtr()[m_outerStart]; } 211 inline Index* innerIndexPtr() 212 { return m_matrix.const_cast_derived().innerIndexPtr() + m_matrix.outerIndexPtr()[m_outerStart]; } 213 214 inline const Index* outerIndexPtr() const 215 { return m_matrix.outerIndexPtr() + m_outerStart; } 216 inline Index* outerIndexPtr() 217 { return m_matrix.const_cast_derived().outerIndexPtr() + m_outerStart; } 218 219 Index nonZeros() const 220 { 221 if(m_matrix.isCompressed()) 222 return std::size_t(m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()]) 223 - std::size_t(m_matrix.outerIndexPtr()[m_outerStart]); 224 else if(m_outerSize.value()==0) 225 return 0; 226 else 227 return Map<const Matrix<Index,OuterSize,1> >(m_matrix.innerNonZeroPtr()+m_outerStart, m_outerSize.value()).sum(); 228 } 229 230 const Scalar& lastCoeff() const 231 { 232 EIGEN_STATIC_ASSERT_VECTOR_ONLY(BlockImpl); 233 eigen_assert(nonZeros()>0); 234 if(m_matrix.isCompressed()) 235 return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart+1]-1]; 236 else 237 return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart]+m_matrix.innerNonZeroPtr()[m_outerStart]-1]; 238 } 239 240 EIGEN_STRONG_INLINE Index rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); } 241 EIGEN_STRONG_INLINE Index cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); } 242 243 protected: 244 245 typename SparseMatrixType::Nested m_matrix; 246 Index m_outerStart; 247 const internal::variable_if_dynamic<Index, OuterSize> m_outerSize; 248 249}; 250 251 252template<typename _Scalar, int _Options, typename _Index, int BlockRows, int BlockCols> 253class BlockImpl<const SparseMatrix<_Scalar, _Options, _Index>,BlockRows,BlockCols,true,Sparse> 254 : public SparseMatrixBase<Block<const SparseMatrix<_Scalar, _Options, _Index>,BlockRows,BlockCols,true> > 255{ 256 typedef SparseMatrix<_Scalar, _Options, _Index> SparseMatrixType; 257 typedef typename internal::remove_all<typename SparseMatrixType::Nested>::type _MatrixTypeNested; 258 typedef Block<const SparseMatrixType, BlockRows, BlockCols, true> BlockType; 259public: 260 enum { IsRowMajor = internal::traits<BlockType>::IsRowMajor }; 261 EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType) 262protected: 263 enum { OuterSize = IsRowMajor ? BlockRows : BlockCols }; 264public: 265 266 class InnerIterator: public SparseMatrixType::InnerIterator 267 { 268 public: 269 inline InnerIterator(const BlockType& xpr, Index outer) 270 : SparseMatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) 271 {} 272 inline Index row() const { return IsRowMajor ? m_outer : this->index(); } 273 inline Index col() const { return IsRowMajor ? this->index() : m_outer; } 274 protected: 275 Index m_outer; 276 }; 277 class ReverseInnerIterator: public SparseMatrixType::ReverseInnerIterator 278 { 279 public: 280 inline ReverseInnerIterator(const BlockType& xpr, Index outer) 281 : SparseMatrixType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) 282 {} 283 inline Index row() const { return IsRowMajor ? m_outer : this->index(); } 284 inline Index col() const { return IsRowMajor ? this->index() : m_outer; } 285 protected: 286 Index m_outer; 287 }; 288 289 inline BlockImpl(const SparseMatrixType& xpr, int i) 290 : m_matrix(xpr), m_outerStart(i), m_outerSize(OuterSize) 291 {} 292 293 inline BlockImpl(const SparseMatrixType& xpr, int startRow, int startCol, int blockRows, int blockCols) 294 : m_matrix(xpr), m_outerStart(IsRowMajor ? startRow : startCol), m_outerSize(IsRowMajor ? blockRows : blockCols) 295 {} 296 297 inline const Scalar* valuePtr() const 298 { return m_matrix.valuePtr() + m_matrix.outerIndexPtr()[m_outerStart]; } 299 300 inline const Index* innerIndexPtr() const 301 { return m_matrix.innerIndexPtr() + m_matrix.outerIndexPtr()[m_outerStart]; } 302 303 inline const Index* outerIndexPtr() const 304 { return m_matrix.outerIndexPtr() + m_outerStart; } 305 306 Index nonZeros() const 307 { 308 if(m_matrix.isCompressed()) 309 return std::size_t(m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()]) 310 - std::size_t(m_matrix.outerIndexPtr()[m_outerStart]); 311 else if(m_outerSize.value()==0) 312 return 0; 313 else 314 return Map<const Matrix<Index,OuterSize,1> >(m_matrix.innerNonZeroPtr()+m_outerStart, m_outerSize.value()).sum(); 315 } 316 317 const Scalar& lastCoeff() const 318 { 319 EIGEN_STATIC_ASSERT_VECTOR_ONLY(BlockImpl); 320 eigen_assert(nonZeros()>0); 321 if(m_matrix.isCompressed()) 322 return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart+1]-1]; 323 else 324 return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart]+m_matrix.innerNonZeroPtr()[m_outerStart]-1]; 325 } 326 327 EIGEN_STRONG_INLINE Index rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); } 328 EIGEN_STRONG_INLINE Index cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); } 329 330 protected: 331 332 typename SparseMatrixType::Nested m_matrix; 333 Index m_outerStart; 334 const internal::variable_if_dynamic<Index, OuterSize> m_outerSize; 335 336}; 337 338//---------- 339 340/** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this 341 * is col-major (resp. row-major). 342 */ 343template<typename Derived> 344typename SparseMatrixBase<Derived>::InnerVectorReturnType SparseMatrixBase<Derived>::innerVector(Index outer) 345{ return InnerVectorReturnType(derived(), outer); } 346 347/** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this 348 * is col-major (resp. row-major). Read-only. 349 */ 350template<typename Derived> 351const typename SparseMatrixBase<Derived>::ConstInnerVectorReturnType SparseMatrixBase<Derived>::innerVector(Index outer) const 352{ return ConstInnerVectorReturnType(derived(), outer); } 353 354/** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this 355 * is col-major (resp. row-major). 356 */ 357template<typename Derived> 358Block<Derived,Dynamic,Dynamic,true> SparseMatrixBase<Derived>::innerVectors(Index outerStart, Index outerSize) 359{ 360 return Block<Derived,Dynamic,Dynamic,true>(derived(), 361 IsRowMajor ? outerStart : 0, IsRowMajor ? 0 : outerStart, 362 IsRowMajor ? outerSize : rows(), IsRowMajor ? cols() : outerSize); 363 364} 365 366/** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this 367 * is col-major (resp. row-major). Read-only. 368 */ 369template<typename Derived> 370const Block<const Derived,Dynamic,Dynamic,true> SparseMatrixBase<Derived>::innerVectors(Index outerStart, Index outerSize) const 371{ 372 return Block<const Derived,Dynamic,Dynamic,true>(derived(), 373 IsRowMajor ? outerStart : 0, IsRowMajor ? 0 : outerStart, 374 IsRowMajor ? outerSize : rows(), IsRowMajor ? cols() : outerSize); 375 376} 377 378/** Generic implementation of sparse Block expression. 379 * Real-only. 380 */ 381template<typename XprType, int BlockRows, int BlockCols, bool InnerPanel> 382class BlockImpl<XprType,BlockRows,BlockCols,InnerPanel,Sparse> 383 : public SparseMatrixBase<Block<XprType,BlockRows,BlockCols,InnerPanel> >, internal::no_assignment_operator 384{ 385 typedef typename internal::remove_all<typename XprType::Nested>::type _MatrixTypeNested; 386 typedef Block<XprType, BlockRows, BlockCols, InnerPanel> BlockType; 387public: 388 enum { IsRowMajor = internal::traits<BlockType>::IsRowMajor }; 389 EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType) 390 391 /** Column or Row constructor 392 */ 393 inline BlockImpl(const XprType& xpr, int i) 394 : m_matrix(xpr), 395 m_startRow( (BlockRows==1) && (BlockCols==XprType::ColsAtCompileTime) ? i : 0), 396 m_startCol( (BlockRows==XprType::RowsAtCompileTime) && (BlockCols==1) ? i : 0), 397 m_blockRows(xpr.rows()), 398 m_blockCols(xpr.cols()) 399 {} 400 401 /** Dynamic-size constructor 402 */ 403 inline BlockImpl(const XprType& xpr, int startRow, int startCol, int blockRows, int blockCols) 404 : m_matrix(xpr), m_startRow(startRow), m_startCol(startCol), m_blockRows(blockRows), m_blockCols(blockCols) 405 {} 406 407 inline int rows() const { return m_blockRows.value(); } 408 inline int cols() const { return m_blockCols.value(); } 409 410 inline Scalar& coeffRef(int row, int col) 411 { 412 return m_matrix.const_cast_derived() 413 .coeffRef(row + m_startRow.value(), col + m_startCol.value()); 414 } 415 416 inline const Scalar coeff(int row, int col) const 417 { 418 return m_matrix.coeff(row + m_startRow.value(), col + m_startCol.value()); 419 } 420 421 inline Scalar& coeffRef(int index) 422 { 423 return m_matrix.const_cast_derived() 424 .coeffRef(m_startRow.value() + (RowsAtCompileTime == 1 ? 0 : index), 425 m_startCol.value() + (RowsAtCompileTime == 1 ? index : 0)); 426 } 427 428 inline const Scalar coeff(int index) const 429 { 430 return m_matrix 431 .coeff(m_startRow.value() + (RowsAtCompileTime == 1 ? 0 : index), 432 m_startCol.value() + (RowsAtCompileTime == 1 ? index : 0)); 433 } 434 435 inline const _MatrixTypeNested& nestedExpression() const { return m_matrix; } 436 437 class InnerIterator : public _MatrixTypeNested::InnerIterator 438 { 439 typedef typename _MatrixTypeNested::InnerIterator Base; 440 const BlockType& m_block; 441 Index m_end; 442 public: 443 444 EIGEN_STRONG_INLINE InnerIterator(const BlockType& block, Index outer) 445 : Base(block.derived().nestedExpression(), outer + (IsRowMajor ? block.m_startRow.value() : block.m_startCol.value())), 446 m_block(block), 447 m_end(IsRowMajor ? block.m_startCol.value()+block.m_blockCols.value() : block.m_startRow.value()+block.m_blockRows.value()) 448 { 449 while( (Base::operator bool()) && (Base::index() < (IsRowMajor ? m_block.m_startCol.value() : m_block.m_startRow.value())) ) 450 Base::operator++(); 451 } 452 453 inline Index index() const { return Base::index() - (IsRowMajor ? m_block.m_startCol.value() : m_block.m_startRow.value()); } 454 inline Index outer() const { return Base::outer() - (IsRowMajor ? m_block.m_startRow.value() : m_block.m_startCol.value()); } 455 inline Index row() const { return Base::row() - m_block.m_startRow.value(); } 456 inline Index col() const { return Base::col() - m_block.m_startCol.value(); } 457 458 inline operator bool() const { return Base::operator bool() && Base::index() < m_end; } 459 }; 460 class ReverseInnerIterator : public _MatrixTypeNested::ReverseInnerIterator 461 { 462 typedef typename _MatrixTypeNested::ReverseInnerIterator Base; 463 const BlockType& m_block; 464 Index m_begin; 465 public: 466 467 EIGEN_STRONG_INLINE ReverseInnerIterator(const BlockType& block, Index outer) 468 : Base(block.derived().nestedExpression(), outer + (IsRowMajor ? block.m_startRow.value() : block.m_startCol.value())), 469 m_block(block), 470 m_begin(IsRowMajor ? block.m_startCol.value() : block.m_startRow.value()) 471 { 472 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())) ) 473 Base::operator--(); 474 } 475 476 inline Index index() const { return Base::index() - (IsRowMajor ? m_block.m_startCol.value() : m_block.m_startRow.value()); } 477 inline Index outer() const { return Base::outer() - (IsRowMajor ? m_block.m_startRow.value() : m_block.m_startCol.value()); } 478 inline Index row() const { return Base::row() - m_block.m_startRow.value(); } 479 inline Index col() const { return Base::col() - m_block.m_startCol.value(); } 480 481 inline operator bool() const { return Base::operator bool() && Base::index() >= m_begin; } 482 }; 483 protected: 484 friend class InnerIterator; 485 friend class ReverseInnerIterator; 486 487 EIGEN_INHERIT_ASSIGNMENT_OPERATORS(BlockImpl) 488 489 typename XprType::Nested m_matrix; 490 const internal::variable_if_dynamic<Index, XprType::RowsAtCompileTime == 1 ? 0 : Dynamic> m_startRow; 491 const internal::variable_if_dynamic<Index, XprType::ColsAtCompileTime == 1 ? 0 : Dynamic> m_startCol; 492 const internal::variable_if_dynamic<Index, RowsAtCompileTime> m_blockRows; 493 const internal::variable_if_dynamic<Index, ColsAtCompileTime> m_blockCols; 494 495}; 496 497} // end namespace Eigen 498 499#endif // EIGEN_SPARSE_BLOCK_H 500