12b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// This file is part of Eigen, a lightweight C++ template library 22b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// for linear algebra. 32b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// 42b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Copyright (C) 2013 Desire Nuentsa <desire.nuentsa_wakam@inria.fr> 52b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Copyright (C) 2013 Gael Guennebaud <gael.guennebaud@inria.fr> 62b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// 72b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// This Source Code Form is subject to the terms of the Mozilla 82b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Public License v. 2.0. If a copy of the MPL was not distributed 92b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#ifndef EIGEN_SPARSEBLOCKMATRIX_H 122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#define EIGEN_SPARSEBLOCKMATRIX_H 132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangnamespace Eigen { 152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang/** \ingroup SparseCore_Module 162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \class BlockSparseMatrix 182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \brief A versatile sparse matrix representation where each element is a block 202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * This class provides routines to manipulate block sparse matrices stored in a 222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * BSR-like representation. There are two main types : 232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 1. All blocks have the same number of rows and columns, called block size 252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * in the following. In this case, if this block size is known at compile time, 262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * it can be given as a template parameter like 272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \code 282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * BlockSparseMatrix<Scalar, 3, ColMajor> bmat(b_rows, b_cols); 292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \endcode 302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * Here, bmat is a b_rows x b_cols block sparse matrix 312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * where each coefficient is a 3x3 dense matrix. 322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * If the block size is fixed but will be given at runtime, 332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \code 342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * BlockSparseMatrix<Scalar, Dynamic, ColMajor> bmat(b_rows, b_cols); 352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * bmat.setBlockSize(block_size); 362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \endcode 372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 2. The second case is for variable-block sparse matrices. 392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * Here each block has its own dimensions. The only restriction is that all the blocks 402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * in a row (resp. a column) should have the same number of rows (resp. of columns). 412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * It is thus required in this case to describe the layout of the matrix by calling 422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * setBlockLayout(rowBlocks, colBlocks). 432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * In any of the previous case, the matrix can be filled by calling setFromTriplets(). 452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * A regular sparse matrix can be converted to a block sparse matrix and vice versa. 462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * It is obviously required to describe the block layout beforehand by calling either 472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * setBlockSize() for fixed-size blocks or setBlockLayout for variable-size blocks. 482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \tparam _Scalar The Scalar type 502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \tparam _BlockAtCompileTime The block layout option. It takes the following values 512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * Dynamic : block size known at runtime 522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * a numeric number : fixed-size block known at compile time 532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename _Scalar, int _BlockAtCompileTime=Dynamic, int _Options=ColMajor, typename _StorageIndex=int> class BlockSparseMatrix; 552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename BlockSparseMatrixT> class BlockSparseMatrixView; 572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangnamespace internal { 592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename _Scalar, int _BlockAtCompileTime, int _Options, typename _Index> 602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangstruct traits<BlockSparseMatrix<_Scalar,_BlockAtCompileTime,_Options, _Index> > 612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef _Scalar Scalar; 632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef _Index Index; 642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Sparse StorageKind; // FIXME Where is it used ?? 652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef MatrixXpr XprKind; 662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang enum { 672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang RowsAtCompileTime = Dynamic, 682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ColsAtCompileTime = Dynamic, 692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang MaxRowsAtCompileTime = Dynamic, 702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang MaxColsAtCompileTime = Dynamic, 712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang BlockSize = _BlockAtCompileTime, 722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Flags = _Options | NestByRefBit | LvalueBit, 732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang CoeffReadCost = NumTraits<Scalar>::ReadCost, 742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang SupportedAccessPatterns = InnerRandomAccessPattern 752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang }; 762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename BlockSparseMatrixT> 782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangstruct traits<BlockSparseMatrixView<BlockSparseMatrixT> > 792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Ref<Matrix<typename BlockSparseMatrixT::Scalar, BlockSparseMatrixT::BlockSize, BlockSparseMatrixT::BlockSize> > Scalar; 812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Ref<Matrix<typename BlockSparseMatrixT::RealScalar, BlockSparseMatrixT::BlockSize, BlockSparseMatrixT::BlockSize> > RealScalar; 822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Function object to sort a triplet list 862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename Iterator, bool IsColMajor> 872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangstruct TripletComp 882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename Iterator::value_type Triplet; 902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang bool operator()(const Triplet& a, const Triplet& b) 912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { if(IsColMajor) 922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return ((a.col() == b.col() && a.row() < b.row()) || (a.col() < b.col())); 932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang else 942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return ((a.row() == b.row() && a.col() < b.col()) || (a.row() < b.row())); 952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} // end namespace internal 982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang/* Proxy to view the block sparse matrix as a regular sparse matrix */ 1012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename BlockSparseMatrixT> 1022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangclass BlockSparseMatrixView : public SparseMatrixBase<BlockSparseMatrixT> 1032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 1042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang public: 1052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Ref<typename BlockSparseMatrixT::BlockScalar> Scalar; 1062b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Ref<typename BlockSparseMatrixT::BlockRealScalar> RealScalar; 1072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename BlockSparseMatrixT::Index Index; 1082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef BlockSparseMatrixT Nested; 1092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang enum { 1102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Flags = BlockSparseMatrixT::Options, 1112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Options = BlockSparseMatrixT::Options, 1122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang RowsAtCompileTime = BlockSparseMatrixT::RowsAtCompileTime, 1132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ColsAtCompileTime = BlockSparseMatrixT::ColsAtCompileTime, 1142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang MaxColsAtCompileTime = BlockSparseMatrixT::MaxColsAtCompileTime, 1152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang MaxRowsAtCompileTime = BlockSparseMatrixT::MaxRowsAtCompileTime 1162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang }; 1172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang public: 1182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang BlockSparseMatrixView(const BlockSparseMatrixT& spblockmat) 1192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : m_spblockmat(spblockmat) 1202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang {} 1212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index outerSize() const 1232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 1242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return (Flags&RowMajorBit) == 1 ? this->rows() : this->cols(); 1252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 1262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index cols() const 1272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 1282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return m_spblockmat.blockCols(); 1292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 1302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index rows() const 1312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 1322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return m_spblockmat.blockRows(); 1332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 1342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Scalar coeff(Index row, Index col) 1352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 1362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return m_spblockmat.coeff(row, col); 1372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 1382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Scalar coeffRef(Index row, Index col) 1392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 1402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return m_spblockmat.coeffRef(row, col); 1412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 1422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Wrapper to iterate over all blocks 1432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang class InnerIterator : public BlockSparseMatrixT::BlockInnerIterator 1442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 1452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang public: 1462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang InnerIterator(const BlockSparseMatrixView& mat, Index outer) 1472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : BlockSparseMatrixT::BlockInnerIterator(mat.m_spblockmat, outer) 1482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang {} 1492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang }; 1512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang protected: 1532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const BlockSparseMatrixT& m_spblockmat; 1542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 1552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Proxy to view a regular vector as a block vector 1572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename BlockSparseMatrixT, typename VectorType> 1582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangclass BlockVectorView 1592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 1602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang public: 1612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang enum { 1622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang BlockSize = BlockSparseMatrixT::BlockSize, 1632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ColsAtCompileTime = VectorType::ColsAtCompileTime, 1642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang RowsAtCompileTime = VectorType::RowsAtCompileTime, 1652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Flags = VectorType::Flags 1662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang }; 1672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Ref<const Matrix<typename BlockSparseMatrixT::Scalar, (RowsAtCompileTime==1)? 1 : BlockSize, (ColsAtCompileTime==1)? 1 : BlockSize> >Scalar; 1682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename BlockSparseMatrixT::Index Index; 1692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang public: 1702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang BlockVectorView(const BlockSparseMatrixT& spblockmat, const VectorType& vec) 1712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : m_spblockmat(spblockmat),m_vec(vec) 1722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { } 1732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index cols() const 1742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 1752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return m_vec.cols(); 1762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 1772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index size() const 1782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 1792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return m_spblockmat.blockRows(); 1802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 1812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Scalar coeff(Index bi) const 1822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 1832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index startRow = m_spblockmat.blockRowsIndex(bi); 1842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index rowSize = m_spblockmat.blockRowsIndex(bi+1) - startRow; 1852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return m_vec.middleRows(startRow, rowSize); 1862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 1872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Scalar coeff(Index bi, Index j) const 1882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 1892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index startRow = m_spblockmat.blockRowsIndex(bi); 1902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index rowSize = m_spblockmat.blockRowsIndex(bi+1) - startRow; 1912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return m_vec.block(startRow, j, rowSize, 1); 1922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 1932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang protected: 1942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const BlockSparseMatrixT& m_spblockmat; 1952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const VectorType& m_vec; 1962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 1972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 1982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename VectorType, typename Index> class BlockVectorReturn; 1992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Proxy to view a regular vector as a block vector 2022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename BlockSparseMatrixT, typename VectorType> 2032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangclass BlockVectorReturn 2042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 2052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang public: 2062b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang enum { 2072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ColsAtCompileTime = VectorType::ColsAtCompileTime, 2082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang RowsAtCompileTime = VectorType::RowsAtCompileTime, 2092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Flags = VectorType::Flags 2102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang }; 2112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Ref<Matrix<typename VectorType::Scalar, RowsAtCompileTime, ColsAtCompileTime> > Scalar; 2122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename BlockSparseMatrixT::Index Index; 2132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang public: 2142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang BlockVectorReturn(const BlockSparseMatrixT& spblockmat, VectorType& vec) 2152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : m_spblockmat(spblockmat),m_vec(vec) 2162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { } 2172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index size() const 2182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 2192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return m_spblockmat.blockRows(); 2202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 2212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Scalar coeffRef(Index bi) 2222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 2232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index startRow = m_spblockmat.blockRowsIndex(bi); 2242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index rowSize = m_spblockmat.blockRowsIndex(bi+1) - startRow; 2252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return m_vec.middleRows(startRow, rowSize); 2262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 2272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Scalar coeffRef(Index bi, Index j) 2282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 2292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index startRow = m_spblockmat.blockRowsIndex(bi); 2302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index rowSize = m_spblockmat.blockRowsIndex(bi+1) - startRow; 2312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return m_vec.block(startRow, j, rowSize, 1); 2322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 2332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang protected: 2352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const BlockSparseMatrixT& m_spblockmat; 2362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VectorType& m_vec; 2372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 2382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Block version of the sparse dense product 2402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename Lhs, typename Rhs> 2412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangclass BlockSparseTimeDenseProduct; 2422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangnamespace internal { 2442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename BlockSparseMatrixT, typename VecType> 2462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangstruct traits<BlockSparseTimeDenseProduct<BlockSparseMatrixT, VecType> > 2472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 2482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Dense StorageKind; 2492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef MatrixXpr XprKind; 2502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename BlockSparseMatrixT::Scalar Scalar; 2512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename BlockSparseMatrixT::Index Index; 2522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang enum { 2532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang RowsAtCompileTime = Dynamic, 2542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ColsAtCompileTime = Dynamic, 2552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang MaxRowsAtCompileTime = Dynamic, 2562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang MaxColsAtCompileTime = Dynamic, 2572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Flags = 0, 2582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang CoeffReadCost = internal::traits<BlockSparseMatrixT>::CoeffReadCost 2592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang }; 2602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 2612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} // end namespace internal 2622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename Lhs, typename Rhs> 2642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangclass BlockSparseTimeDenseProduct 2652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : public ProductBase<BlockSparseTimeDenseProduct<Lhs,Rhs>, Lhs, Rhs> 2662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 2672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang public: 2682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang EIGEN_PRODUCT_PUBLIC_INTERFACE(BlockSparseTimeDenseProduct) 2692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang BlockSparseTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) 2712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang {} 2722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename Dest> void scaleAndAddTo(Dest& dest, const typename Rhs::Scalar& alpha) const 2742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 2752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang BlockVectorReturn<Lhs,Dest> tmpDest(m_lhs, dest); 2762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang internal::sparse_time_dense_product( BlockSparseMatrixView<Lhs>(m_lhs), BlockVectorView<Lhs, Rhs>(m_lhs, m_rhs), tmpDest, alpha); 2772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 2782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang private: 2802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang BlockSparseTimeDenseProduct& operator=(const BlockSparseTimeDenseProduct&); 2812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 2822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename _Scalar, int _BlockAtCompileTime, int _Options, typename _StorageIndex> 2842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangclass BlockSparseMatrix : public SparseMatrixBase<BlockSparseMatrix<_Scalar,_BlockAtCompileTime, _Options,_StorageIndex> > 2852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 2862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang public: 2872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef _Scalar Scalar; 2882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename NumTraits<Scalar>::Real RealScalar; 2892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef _StorageIndex StorageIndex; 2902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename internal::ref_selector<BlockSparseMatrix<_Scalar, _BlockAtCompileTime, _Options, _StorageIndex> >::type Nested; 2912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 2922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang enum { 2932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Options = _Options, 2942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Flags = Options, 2952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang BlockSize=_BlockAtCompileTime, 2962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang RowsAtCompileTime = Dynamic, 2972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ColsAtCompileTime = Dynamic, 2982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang MaxRowsAtCompileTime = Dynamic, 2992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang MaxColsAtCompileTime = Dynamic, 3002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang IsVectorAtCompileTime = 0, 3012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang IsColMajor = Flags&RowMajorBit ? 0 : 1 3022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang }; 3032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Matrix<Scalar, _BlockAtCompileTime, _BlockAtCompileTime,IsColMajor ? ColMajor : RowMajor> BlockScalar; 3042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Matrix<RealScalar, _BlockAtCompileTime, _BlockAtCompileTime,IsColMajor ? ColMajor : RowMajor> BlockRealScalar; 3052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef typename internal::conditional<_BlockAtCompileTime==Dynamic, Scalar, BlockScalar>::type BlockScalarReturnType; 3062b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef BlockSparseMatrix<Scalar, BlockSize, IsColMajor ? ColMajor : RowMajor, StorageIndex> PlainObject; 3072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang public: 3082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Default constructor 3092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang BlockSparseMatrix() 3102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : m_innerBSize(0),m_outerBSize(0),m_innerOffset(0),m_outerOffset(0), 3112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_nonzerosblocks(0),m_values(0),m_blockPtr(0),m_indices(0), 3122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_outerIndex(0),m_blockSize(BlockSize) 3132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { } 3142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** 3172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \brief Construct and resize 3182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 3192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 3202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang BlockSparseMatrix(Index brow, Index bcol) 3212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : m_innerBSize(IsColMajor ? brow : bcol), 3222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_outerBSize(IsColMajor ? bcol : brow), 3232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_innerOffset(0),m_outerOffset(0),m_nonzerosblocks(0), 3242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_values(0),m_blockPtr(0),m_indices(0), 3252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_outerIndex(0),m_blockSize(BlockSize) 3262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { } 3272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** 3292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \brief Copy-constructor 3302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 3312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang BlockSparseMatrix(const BlockSparseMatrix& other) 3322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : m_innerBSize(other.m_innerBSize),m_outerBSize(other.m_outerBSize), 3332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_nonzerosblocks(other.m_nonzerosblocks),m_nonzeros(other.m_nonzeros), 3342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_blockPtr(0),m_blockSize(other.m_blockSize) 3352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 3362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // should we allow copying between variable-size blocks and fixed-size blocks ?? 3372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(m_blockSize == BlockSize && " CAN NOT COPY BETWEEN FIXED-SIZE AND VARIABLE-SIZE BLOCKS"); 3382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::copy(other.m_innerOffset, other.m_innerOffset+m_innerBSize+1, m_innerOffset); 3402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::copy(other.m_outerOffset, other.m_outerOffset+m_outerBSize+1, m_outerOffset); 3412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::copy(other.m_values, other.m_values+m_nonzeros, m_values); 3422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(m_blockSize != Dynamic) 3442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::copy(other.m_blockPtr, other.m_blockPtr+m_nonzerosblocks, m_blockPtr); 3452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::copy(other.m_indices, other.m_indices+m_nonzerosblocks, m_indices); 3472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::copy(other.m_outerIndex, other.m_outerIndex+m_outerBSize, m_outerIndex); 3482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 3492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang friend void swap(BlockSparseMatrix& first, BlockSparseMatrix& second) 3512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 3522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::swap(first.m_innerBSize, second.m_innerBSize); 3532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::swap(first.m_outerBSize, second.m_outerBSize); 3542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::swap(first.m_innerOffset, second.m_innerOffset); 3552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::swap(first.m_outerOffset, second.m_outerOffset); 3562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::swap(first.m_nonzerosblocks, second.m_nonzerosblocks); 3572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::swap(first.m_nonzeros, second.m_nonzeros); 3582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::swap(first.m_values, second.m_values); 3592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::swap(first.m_blockPtr, second.m_blockPtr); 3602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::swap(first.m_indices, second.m_indices); 3612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::swap(first.m_outerIndex, second.m_outerIndex); 3622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::swap(first.m_BlockSize, second.m_blockSize); 3632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 3642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang BlockSparseMatrix& operator=(BlockSparseMatrix other) 3662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 3672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang //Copy-and-swap paradigm ... avoid leaked data if thrown 3682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang swap(*this, other); 3692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return *this; 3702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 3712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Destructor 3732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ~BlockSparseMatrix() 3742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 3752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang delete[] m_outerIndex; 3762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang delete[] m_innerOffset; 3772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang delete[] m_outerOffset; 3782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang delete[] m_indices; 3792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang delete[] m_blockPtr; 3802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang delete[] m_values; 3812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 3822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** 3852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \brief Constructor from a sparse matrix 3862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 3872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 3882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename MatrixType> 3892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline BlockSparseMatrix(const MatrixType& spmat) : m_blockSize(BlockSize) 3902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 3912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang EIGEN_STATIC_ASSERT((m_blockSize != Dynamic), THIS_METHOD_IS_ONLY_FOR_FIXED_SIZE); 3922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang *this = spmat; 3942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 3952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** 3972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \brief Assignment from a sparse matrix with the same storage order 3982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 3992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * Convert from a sparse matrix to block sparse matrix. 4002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \warning Before calling this function, tt is necessary to call 4012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * either setBlockLayout() (matrices with variable-size blocks) 4022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * or setBlockSize() (for fixed-size blocks). 4032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 4042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename MatrixType> 4052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline BlockSparseMatrix& operator=(const MatrixType& spmat) 4062b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 4072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert((m_innerBSize != 0 && m_outerBSize != 0) 4082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang && "Trying to assign to a zero-size matrix, call resize() first"); 4092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(((MatrixType::Options&RowMajorBit) != IsColMajor) && "Wrong storage order"); 4102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef SparseMatrix<bool,MatrixType::Options,typename MatrixType::Index> MatrixPatternType; 4112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang MatrixPatternType blockPattern(blockRows(), blockCols()); 4122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_nonzeros = 0; 4132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 4142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // First, compute the number of nonzero blocks and their locations 4152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(StorageIndex bj = 0; bj < m_outerBSize; ++bj) 4162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 4172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Browse each outer block and compute the structure 4182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::vector<bool> nzblocksFlag(m_innerBSize,false); // Record the existing blocks 4192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang blockPattern.startVec(bj); 4202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(StorageIndex j = blockOuterIndex(bj); j < blockOuterIndex(bj+1); ++j) 4212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 4222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typename MatrixType::InnerIterator it_spmat(spmat, j); 4232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(; it_spmat; ++it_spmat) 4242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 4252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex bi = innerToBlock(it_spmat.index()); // Index of the current nonzero block 4262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(!nzblocksFlag[bi]) 4272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 4282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Save the index of this nonzero block 4292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang nzblocksFlag[bi] = true; 4302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang blockPattern.insertBackByOuterInnerUnordered(bj, bi) = true; 4312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Compute the total number of nonzeros (including explicit zeros in blocks) 4322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_nonzeros += blockOuterSize(bj) * blockInnerSize(bi); 4332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 4342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 4352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } // end current outer block 4362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 4372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang blockPattern.finalize(); 4382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 4392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Allocate the internal arrays 4402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang setBlockStructure(blockPattern); 4412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 4422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(StorageIndex nz = 0; nz < m_nonzeros; ++nz) m_values[nz] = Scalar(0); 4432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(StorageIndex bj = 0; bj < m_outerBSize; ++bj) 4442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 4452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Now copy the values 4462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(StorageIndex j = blockOuterIndex(bj); j < blockOuterIndex(bj+1); ++j) 4472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 4482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Browse the outer block column by column (for column-major matrices) 4492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typename MatrixType::InnerIterator it_spmat(spmat, j); 4502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(; it_spmat; ++it_spmat) 4512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 4522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex idx = 0; // Position of this block in the column block 4532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex bi = innerToBlock(it_spmat.index()); // Index of the current nonzero block 4542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Go to the inner block where this element belongs to 4552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang while(bi > m_indices[m_outerIndex[bj]+idx]) ++idx; // Not expensive for ordered blocks 4562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex idxVal;// Get the right position in the array of values for this element 4572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(m_blockSize == Dynamic) 4582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 4592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Offset from all blocks before ... 4602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang idxVal = m_blockPtr[m_outerIndex[bj]+idx]; 4612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // ... and offset inside the block 4622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang idxVal += (j - blockOuterIndex(bj)) * blockOuterSize(bj) + it_spmat.index() - m_innerOffset[bi]; 4632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 4642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang else 4652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 4662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // All blocks before 4672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang idxVal = (m_outerIndex[bj] + idx) * m_blockSize * m_blockSize; 4682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // inside the block 4692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang idxVal += (j - blockOuterIndex(bj)) * m_blockSize + (it_spmat.index()%m_blockSize); 4702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 4712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Insert the value 4722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_values[idxVal] = it_spmat.value(); 4732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } // end of this column 4742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } // end of this block 4752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } // end of this outer block 4762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 4772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return *this; 4782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 4792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 4802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** 4812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \brief Set the nonzero block pattern of the matrix 4822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 4832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * Given a sparse matrix describing the nonzero block pattern, 4842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * this function prepares the internal pointers for values. 4852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * After calling this function, any *nonzero* block (bi, bj) can be set 4862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * with a simple call to coeffRef(bi,bj). 4872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 4882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 4892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \warning Before calling this function, tt is necessary to call 4902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * either setBlockLayout() (matrices with variable-size blocks) 4912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * or setBlockSize() (for fixed-size blocks). 4922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 4932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \param blockPattern Sparse matrix of boolean elements describing the block structure 4942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 4952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \sa setBlockLayout() \sa setBlockSize() 4962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 4972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename MatrixType> 4982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang void setBlockStructure(const MatrixType& blockPattern) 4992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 5002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang resize(blockPattern.rows(), blockPattern.cols()); 5012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang reserve(blockPattern.nonZeros()); 5022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 5032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Browse the block pattern and set up the various pointers 5042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_outerIndex[0] = 0; 5052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(m_blockSize == Dynamic) m_blockPtr[0] = 0; 5062b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(StorageIndex nz = 0; nz < m_nonzeros; ++nz) m_values[nz] = Scalar(0); 5072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(StorageIndex bj = 0; bj < m_outerBSize; ++bj) 5082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 5092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang //Browse each outer block 5102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 5112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang //First, copy and save the indices of nonzero blocks 5122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang //FIXME : find a way to avoid this ... 5132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::vector<int> nzBlockIdx; 5142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typename MatrixType::InnerIterator it(blockPattern, bj); 5152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(; it; ++it) 5162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 5172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang nzBlockIdx.push_back(it.index()); 5182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 5192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::sort(nzBlockIdx.begin(), nzBlockIdx.end()); 5202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 5212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Now, fill block indices and (eventually) pointers to blocks 5222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(StorageIndex idx = 0; idx < nzBlockIdx.size(); ++idx) 5232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 5242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex offset = m_outerIndex[bj]+idx; // offset in m_indices 5252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_indices[offset] = nzBlockIdx[idx]; 5262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(m_blockSize == Dynamic) 5272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_blockPtr[offset] = m_blockPtr[offset-1] + blockInnerSize(nzBlockIdx[idx]) * blockOuterSize(bj); 5282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // There is no blockPtr for fixed-size blocks... not needed !??? 5292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 5302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Save the pointer to the next outer block 5312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_outerIndex[bj+1] = m_outerIndex[bj] + nzBlockIdx.size(); 5322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 5332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 5342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 5352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** 5362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \brief Set the number of rows and columns blocks 5372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 5382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline void resize(Index brow, Index bcol) 5392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 5402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_innerBSize = IsColMajor ? brow : bcol; 5412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_outerBSize = IsColMajor ? bcol : brow; 5422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 5432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 5442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** 5452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \brief set the block size at runtime for fixed-size block layout 5462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 5472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * Call this only for fixed-size blocks 5482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 5492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline void setBlockSize(Index blockSize) 5502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 5512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_blockSize = blockSize; 5522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 5532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 5542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** 5552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \brief Set the row and column block layouts, 5562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 5572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * This function set the size of each row and column block. 5582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * So this function should be used only for blocks with variable size. 5592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \param rowBlocks : Number of rows per row block 5602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \param colBlocks : Number of columns per column block 5612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \sa resize(), setBlockSize() 5622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 5632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline void setBlockLayout(const VectorXi& rowBlocks, const VectorXi& colBlocks) 5642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 5652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const VectorXi& innerBlocks = IsColMajor ? rowBlocks : colBlocks; 5662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const VectorXi& outerBlocks = IsColMajor ? colBlocks : rowBlocks; 5672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(m_innerBSize == innerBlocks.size() && "CHECK THE NUMBER OF ROW OR COLUMN BLOCKS"); 5682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(m_outerBSize == outerBlocks.size() && "CHECK THE NUMBER OF ROW OR COLUMN BLOCKS"); 5692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_outerBSize = outerBlocks.size(); 5702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // starting index of blocks... cumulative sums 5712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_innerOffset = new StorageIndex[m_innerBSize+1]; 5722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_outerOffset = new StorageIndex[m_outerBSize+1]; 5732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_innerOffset[0] = 0; 5742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_outerOffset[0] = 0; 5752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::partial_sum(&innerBlocks[0], &innerBlocks[m_innerBSize-1]+1, &m_innerOffset[1]); 5762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::partial_sum(&outerBlocks[0], &outerBlocks[m_outerBSize-1]+1, &m_outerOffset[1]); 5772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 5782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Compute the total number of nonzeros 5792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_nonzeros = 0; 5802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(StorageIndex bj = 0; bj < m_outerBSize; ++bj) 5812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(StorageIndex bi = 0; bi < m_innerBSize; ++bi) 5822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_nonzeros += outerBlocks[bj] * innerBlocks[bi]; 5832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 5842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 5852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 5862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** 5872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \brief Allocate the internal array of pointers to blocks and their inner indices 5882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 5892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \note For fixed-size blocks, call setBlockSize() to set the block. 5902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * And For variable-size blocks, call setBlockLayout() before using this function 5912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 5922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \param nonzerosblocks Number of nonzero blocks. The total number of nonzeros is 5932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * is computed in setBlockLayout() for variable-size blocks 5942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \sa setBlockSize() 5952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 5962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline void reserve(const Index nonzerosblocks) 5972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 5982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert((m_innerBSize != 0 && m_outerBSize != 0) && 5992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang "TRYING TO RESERVE ZERO-SIZE MATRICES, CALL resize() first"); 6002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 6012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang //FIXME Should free if already allocated 6022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_outerIndex = new StorageIndex[m_outerBSize+1]; 6032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 6042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_nonzerosblocks = nonzerosblocks; 6052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(m_blockSize != Dynamic) 6062b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 6072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_nonzeros = nonzerosblocks * (m_blockSize * m_blockSize); 6082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_blockPtr = 0; 6092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 6102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang else 6112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 6122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // m_nonzeros is already computed in setBlockLayout() 6132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_blockPtr = new StorageIndex[m_nonzerosblocks+1]; 6142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 6152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_indices = new StorageIndex[m_nonzerosblocks+1]; 6162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_values = new Scalar[m_nonzeros]; 6172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 6182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 6192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 6202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** 6212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \brief Fill values in a matrix from a triplet list. 6222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 6232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * Each triplet item has a block stored in an Eigen dense matrix. 6242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * The InputIterator class should provide the functions row(), col() and value() 6252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 6262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \note For fixed-size blocks, call setBlockSize() before this function. 6272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * 6282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * FIXME Do not accept duplicates 6292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 6302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename InputIterator> 6312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang void setFromTriplets(const InputIterator& begin, const InputIterator& end) 6322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 6332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert((m_innerBSize!=0 && m_outerBSize !=0) && "ZERO BLOCKS, PLEASE CALL resize() before"); 6342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 6352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /* First, sort the triplet list 6362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * FIXME This can be unnecessarily expensive since only the inner indices have to be sorted 6372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * The best approach is like in SparseMatrix::setFromTriplets() 6382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 6392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang internal::TripletComp<InputIterator, IsColMajor> tripletcomp; 6402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang std::sort(begin, end, tripletcomp); 6412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 6422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /* Count the number of rows and column blocks, 6432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * and the number of nonzero blocks per outer dimension 6442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 6452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VectorXi rowBlocks(m_innerBSize); // Size of each block row 6462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VectorXi colBlocks(m_outerBSize); // Size of each block column 6472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang rowBlocks.setZero(); colBlocks.setZero(); 6482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VectorXi nzblock_outer(m_outerBSize); // Number of nz blocks per outer vector 6492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VectorXi nz_outer(m_outerBSize); // Number of nz per outer vector...for variable-size blocks 6502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang nzblock_outer.setZero(); 6512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang nz_outer.setZero(); 6522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(InputIterator it(begin); it !=end; ++it) 6532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 6542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(it->row() >= 0 && it->row() < this->blockRows() && it->col() >= 0 && it->col() < this->blockCols()); 6552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert((it->value().rows() == it->value().cols() && (it->value().rows() == m_blockSize)) 6562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang || (m_blockSize == Dynamic)); 6572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 6582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(m_blockSize == Dynamic) 6592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 6602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert((rowBlocks[it->row()] == 0 || rowBlocks[it->row()] == it->value().rows()) && 6612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang "NON CORRESPONDING SIZES FOR ROW BLOCKS"); 6622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert((colBlocks[it->col()] == 0 || colBlocks[it->col()] == it->value().cols()) && 6632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang "NON CORRESPONDING SIZES FOR COLUMN BLOCKS"); 6642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang rowBlocks[it->row()] =it->value().rows(); 6652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang colBlocks[it->col()] = it->value().cols(); 6662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 6672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang nz_outer(IsColMajor ? it->col() : it->row()) += it->value().rows() * it->value().cols(); 6682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang nzblock_outer(IsColMajor ? it->col() : it->row())++; 6692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 6702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Allocate member arrays 6712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(m_blockSize == Dynamic) setBlockLayout(rowBlocks, colBlocks); 6722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex nzblocks = nzblock_outer.sum(); 6732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang reserve(nzblocks); 6742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 6752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Temporary markers 6762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VectorXi block_id(m_outerBSize); // To be used as a block marker during insertion 6772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 6782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Setup outer index pointers and markers 6792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_outerIndex[0] = 0; 6802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (m_blockSize == Dynamic) m_blockPtr[0] = 0; 6812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(StorageIndex bj = 0; bj < m_outerBSize; ++bj) 6822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 6832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_outerIndex[bj+1] = m_outerIndex[bj] + nzblock_outer(bj); 6842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang block_id(bj) = m_outerIndex[bj]; 6852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(m_blockSize==Dynamic) 6862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 6872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_blockPtr[m_outerIndex[bj+1]] = m_blockPtr[m_outerIndex[bj]] + nz_outer(bj); 6882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 6892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 6902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 6912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Fill the matrix 6922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(InputIterator it(begin); it!=end; ++it) 6932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 6942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex outer = IsColMajor ? it->col() : it->row(); 6952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex inner = IsColMajor ? it->row() : it->col(); 6962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_indices[block_id(outer)] = inner; 6972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex block_size = it->value().rows()*it->value().cols(); 6982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex nz_marker = blockPtr(block_id[outer]); 6992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang memcpy(&(m_values[nz_marker]), it->value().data(), block_size * sizeof(Scalar)); 7002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(m_blockSize == Dynamic) 7012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 7022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_blockPtr[block_id(outer)+1] = m_blockPtr[block_id(outer)] + block_size; 7032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 7042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang block_id(outer)++; 7052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 7062b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 7072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // An alternative when the outer indices are sorted...no need to use an array of markers 7082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// for(Index bcol = 0; bcol < m_outerBSize; ++bcol) 7092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// { 7102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Index id = 0, id_nz = 0, id_nzblock = 0; 7112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// for(InputIterator it(begin); it!=end; ++it) 7122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// { 7132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// while (id<bcol) // one pass should do the job unless there are empty columns 7142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// { 7152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// id++; 7162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// m_outerIndex[id+1]=m_outerIndex[id]; 7172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// } 7182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// m_outerIndex[id+1] += 1; 7192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// m_indices[id_nzblock]=brow; 7202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Index block_size = it->value().rows()*it->value().cols(); 7212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// m_blockPtr[id_nzblock+1] = m_blockPtr[id_nzblock] + block_size; 7222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// id_nzblock++; 7232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// memcpy(&(m_values[id_nz]),it->value().data(), block_size*sizeof(Scalar)); 7242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// id_nz += block_size; 7252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// } 7262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// while(id < m_outerBSize-1) // Empty columns at the end 7272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// { 7282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// id++; 7292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// m_outerIndex[id+1]=m_outerIndex[id]; 7302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// } 7312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// } 7322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 7332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 7342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 7352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** 7362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \returns the number of rows 7372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 7382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index rows() const 7392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 7402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// return blockRows(); 7412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return (IsColMajor ? innerSize() : outerSize()); 7422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 7432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 7442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** 7452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \returns the number of cols 7462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 7472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index cols() const 7482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 7492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// return blockCols(); 7502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return (IsColMajor ? outerSize() : innerSize()); 7512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 7522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 7532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index innerSize() const 7542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 7552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(m_blockSize == Dynamic) return m_innerOffset[m_innerBSize]; 7562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang else return (m_innerBSize * m_blockSize) ; 7572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 7582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 7592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index outerSize() const 7602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 7612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(m_blockSize == Dynamic) return m_outerOffset[m_outerBSize]; 7622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang else return (m_outerBSize * m_blockSize) ; 7632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 7642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** \returns the number of rows grouped by blocks */ 7652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index blockRows() const 7662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 7672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return (IsColMajor ? m_innerBSize : m_outerBSize); 7682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 7692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** \returns the number of columns grouped by blocks */ 7702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index blockCols() const 7712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 7722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return (IsColMajor ? m_outerBSize : m_innerBSize); 7732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 7742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 7752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index outerBlocks() const { return m_outerBSize; } 7762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index innerBlocks() const { return m_innerBSize; } 7772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 7782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** \returns the block index where outer belongs to */ 7792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index outerToBlock(Index outer) const 7802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 7812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(outer < outerSize() && "OUTER INDEX OUT OF BOUNDS"); 7822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 7832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(m_blockSize != Dynamic) 7842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return (outer / m_blockSize); // Integer division 7852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 7862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex b_outer = 0; 7872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang while(m_outerOffset[b_outer] <= outer) ++b_outer; 7882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return b_outer - 1; 7892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 7902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** \returns the block index where inner belongs to */ 7912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index innerToBlock(Index inner) const 7922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 7932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(inner < innerSize() && "OUTER INDEX OUT OF BOUNDS"); 7942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 7952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(m_blockSize != Dynamic) 7962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return (inner / m_blockSize); // Integer division 7972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 7982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex b_inner = 0; 7992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang while(m_innerOffset[b_inner] <= inner) ++b_inner; 8002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return b_inner - 1; 8012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 8022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 8032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** 8042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang *\returns a reference to the (i,j) block as an Eigen Dense Matrix 8052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 8062b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Ref<BlockScalar> coeffRef(Index brow, Index bcol) 8072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 8082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(brow < blockRows() && "BLOCK ROW INDEX OUT OF BOUNDS"); 8092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(bcol < blockCols() && "BLOCK nzblocksFlagCOLUMN OUT OF BOUNDS"); 8102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 8112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex rsize = IsColMajor ? blockInnerSize(brow): blockOuterSize(bcol); 8122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex csize = IsColMajor ? blockOuterSize(bcol) : blockInnerSize(brow); 8132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex inner = IsColMajor ? brow : bcol; 8142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex outer = IsColMajor ? bcol : brow; 8152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex offset = m_outerIndex[outer]; 8162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang while(offset < m_outerIndex[outer+1] && m_indices[offset] != inner) 8172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang offset++; 8182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(m_indices[offset] == inner) 8192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 8202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return Map<BlockScalar>(&(m_values[blockPtr(offset)]), rsize, csize); 8212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 8222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang else 8232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 8242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang //FIXME the block does not exist, Insert it !!!!!!!!! 8252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert("DYNAMIC INSERTION IS NOT YET SUPPORTED"); 8262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 8272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 8282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 8292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** 8302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \returns the value of the (i,j) block as an Eigen Dense Matrix 8312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 8322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Map<const BlockScalar> coeff(Index brow, Index bcol) const 8332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 8342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(brow < blockRows() && "BLOCK ROW INDEX OUT OF BOUNDS"); 8352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(bcol < blockCols() && "BLOCK COLUMN OUT OF BOUNDS"); 8362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 8372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex rsize = IsColMajor ? blockInnerSize(brow): blockOuterSize(bcol); 8382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex csize = IsColMajor ? blockOuterSize(bcol) : blockInnerSize(brow); 8392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex inner = IsColMajor ? brow : bcol; 8402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex outer = IsColMajor ? bcol : brow; 8412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex offset = m_outerIndex[outer]; 8422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang while(offset < m_outerIndex[outer+1] && m_indices[offset] != inner) offset++; 8432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(m_indices[offset] == inner) 8442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 8452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return Map<const BlockScalar> (&(m_values[blockPtr(offset)]), rsize, csize); 8462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 8472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang else 8482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// return BlockScalar::Zero(rsize, csize); 8492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert("NOT YET SUPPORTED"); 8502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 8512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 8522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Block Matrix times vector product 8532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename VecType> 8542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang BlockSparseTimeDenseProduct<BlockSparseMatrix, VecType> operator*(const VecType& lhs) const 8552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 8562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return BlockSparseTimeDenseProduct<BlockSparseMatrix, VecType>(*this, lhs); 8572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 8582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 8592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** \returns the number of nonzero blocks */ 8602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index nonZerosBlocks() const { return m_nonzerosblocks; } 8612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** \returns the total number of nonzero elements, including eventual explicit zeros in blocks */ 8622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index nonZeros() const { return m_nonzeros; } 8632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 8642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline BlockScalarReturnType *valuePtr() {return static_cast<BlockScalarReturnType *>(m_values);} 8652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// inline Scalar *valuePtr(){ return m_values; } 8662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline StorageIndex *innerIndexPtr() {return m_indices; } 8672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline const StorageIndex *innerIndexPtr() const {return m_indices; } 8682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline StorageIndex *outerIndexPtr() {return m_outerIndex; } 8692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline const StorageIndex* outerIndexPtr() const {return m_outerIndex; } 8702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 8712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** \brief for compatibility purposes with the SparseMatrix class */ 8722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline bool isCompressed() const {return true;} 8732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** 8742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \returns the starting index of the bi row block 8752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 8762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index blockRowsIndex(Index bi) const 8772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 8782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return IsColMajor ? blockInnerIndex(bi) : blockOuterIndex(bi); 8792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 8802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 8812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** 8822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \returns the starting index of the bj col block 8832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 8842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index blockColsIndex(Index bj) const 8852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 8862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return IsColMajor ? blockOuterIndex(bj) : blockInnerIndex(bj); 8872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 8882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 8892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index blockOuterIndex(Index bj) const 8902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 8912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return (m_blockSize == Dynamic) ? m_outerOffset[bj] : (bj * m_blockSize); 8922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 8932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index blockInnerIndex(Index bi) const 8942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 8952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return (m_blockSize == Dynamic) ? m_innerOffset[bi] : (bi * m_blockSize); 8962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 8972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 8982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Not needed ??? 8992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index blockInnerSize(Index bi) const 9002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 9012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return (m_blockSize == Dynamic) ? (m_innerOffset[bi+1] - m_innerOffset[bi]) : m_blockSize; 9022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 9032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index blockOuterSize(Index bj) const 9042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 9052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return (m_blockSize == Dynamic) ? (m_outerOffset[bj+1]- m_outerOffset[bj]) : m_blockSize; 9062b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 9072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 9082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** 9092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \brief Browse the matrix by outer index 9102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 9112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang class InnerIterator; // Browse column by column 9122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 9132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** 9142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \brief Browse the matrix by block outer index 9152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 9162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang class BlockInnerIterator; // Browse block by block 9172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 9182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang friend std::ostream & operator << (std::ostream & s, const BlockSparseMatrix& m) 9192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 9202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for (StorageIndex j = 0; j < m.outerBlocks(); ++j) 9212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 9222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang BlockInnerIterator itb(m, j); 9232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for(; itb; ++itb) 9242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 9252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang s << "("<<itb.row() << ", " << itb.col() << ")\n"; 9262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang s << itb.value() <<"\n"; 9272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 9282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 9292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang s << std::endl; 9302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return s; 9312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 9322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 9332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** 9342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \returns the starting position of the block <id> in the array of values 9352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang */ 9362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index blockPtr(Index id) const 9372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 9382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(m_blockSize == Dynamic) return m_blockPtr[id]; 9392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang else return id * m_blockSize * m_blockSize; 9402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang //return blockDynIdx(id, typename internal::conditional<(BlockSize==Dynamic), internal::true_type, internal::false_type>::type()); 9412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 9422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 9432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 9442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang protected: 9452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// inline Index blockDynIdx(Index id, internal::true_type) const 9462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// { 9472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// return m_blockPtr[id]; 9482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// } 9492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// inline Index blockDynIdx(Index id, internal::false_type) const 9502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// { 9512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// return id * BlockSize * BlockSize; 9522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// } 9532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 9542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // To be implemented 9552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Insert a block at a particular location... need to make a room for that 9562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Map<BlockScalar> insert(Index brow, Index bcol); 9572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 9582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index m_innerBSize; // Number of block rows 9592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index m_outerBSize; // Number of block columns 9602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex *m_innerOffset; // Starting index of each inner block (size m_innerBSize+1) 9612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex *m_outerOffset; // Starting index of each outer block (size m_outerBSize+1) 9622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index m_nonzerosblocks; // Total nonzeros blocks (lower than m_innerBSize x m_outerBSize) 9632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index m_nonzeros; // Total nonzeros elements 9642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Scalar *m_values; //Values stored block column after block column (size m_nonzeros) 9652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex *m_blockPtr; // Pointer to the beginning of each block in m_values, size m_nonzeroblocks ... null for fixed-size blocks 9662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex *m_indices; //Inner block indices, size m_nonzerosblocks ... OK 9672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang StorageIndex *m_outerIndex; // Starting pointer of each block column in m_indices (size m_outerBSize)... OK 9682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index m_blockSize; // Size of a block for fixed-size blocks, otherwise -1 9692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 9702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 9712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename _Scalar, int _BlockAtCompileTime, int _Options, typename _StorageIndex> 9722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangclass BlockSparseMatrix<_Scalar, _BlockAtCompileTime, _Options, _StorageIndex>::BlockInnerIterator 9732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 9742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang public: 9752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 9762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang enum{ 9772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Flags = _Options 9782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang }; 9792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 9802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang BlockInnerIterator(const BlockSparseMatrix& mat, const Index outer) 9812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : m_mat(mat),m_outer(outer), 9822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_id(mat.m_outerIndex[outer]), 9832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_end(mat.m_outerIndex[outer+1]) 9842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 9852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 9862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 9872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline BlockInnerIterator& operator++() {m_id++; return *this; } 9882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 9892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline const Map<const BlockScalar> value() const 9902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 9912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return Map<const BlockScalar>(&(m_mat.m_values[m_mat.blockPtr(m_id)]), 9922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang rows(),cols()); 9932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 9942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Map<BlockScalar> valueRef() 9952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 9962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return Map<BlockScalar>(&(m_mat.m_values[m_mat.blockPtr(m_id)]), 9972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang rows(),cols()); 9982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 9992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Block inner index 10002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index index() const {return m_mat.m_indices[m_id]; } 10012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index outer() const { return m_outer; } 10022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // block row index 10032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index row() const {return index(); } 10042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // block column index 10052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index col() const {return outer(); } 10062b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // FIXME Number of rows in the current block 10072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 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; } 10082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Number of columns in the current block ... 10092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 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;} 10102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline operator bool() const { return (m_id < m_end); } 10112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 10122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang protected: 10132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const BlockSparseMatrix<_Scalar, _BlockAtCompileTime, _Options, StorageIndex>& m_mat; 10142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const Index m_outer; 10152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index m_id; 10162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index m_end; 10172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 10182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 10192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename _Scalar, int _BlockAtCompileTime, int _Options, typename _StorageIndex> 10202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangclass BlockSparseMatrix<_Scalar, _BlockAtCompileTime, _Options, _StorageIndex>::InnerIterator 10212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 10222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang public: 10232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang InnerIterator(const BlockSparseMatrix& mat, Index outer) 10242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : m_mat(mat),m_outerB(mat.outerToBlock(outer)),m_outer(outer), 10252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang itb(mat, mat.outerToBlock(outer)), 10262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_offset(outer - mat.blockOuterIndex(m_outerB)) 10272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 10282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (itb) 10292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 10302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_id = m_mat.blockInnerIndex(itb.index()); 10312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_start = m_id; 10322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_end = m_mat.blockInnerIndex(itb.index()+1); 10332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 10342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 10352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline InnerIterator& operator++() 10362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 10372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_id++; 10382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (m_id >= m_end) 10392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 10402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ++itb; 10412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (itb) 10422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 10432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_id = m_mat.blockInnerIndex(itb.index()); 10442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_start = m_id; 10452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_end = m_mat.blockInnerIndex(itb.index()+1); 10462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 10472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 10482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return *this; 10492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 10502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline const Scalar& value() const 10512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 10522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return itb.value().coeff(m_id - m_start, m_offset); 10532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 10542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Scalar& valueRef() 10552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 10562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return itb.valueRef().coeff(m_id - m_start, m_offset); 10572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 10582b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index index() const { return m_id; } 10592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index outer() const {return m_outer; } 10602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index col() const {return outer(); } 10612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline Index row() const { return index();} 10622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang inline operator bool() const 10632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 10642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return itb; 10652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 10662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang protected: 10672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const BlockSparseMatrix& m_mat; 10682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const Index m_outer; 10692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const Index m_outerB; 10702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang BlockInnerIterator itb; // Iterator through the blocks 10712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const Index m_offset; // Position of this column in the block 10722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index m_start; // starting inner index of this block 10732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index m_id; // current inner index in the block 10742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index m_end; // starting inner index of the next block 10752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 10762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang}; 10772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} // end namespace Eigen 10782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 10792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#endif // EIGEN_SPARSEBLOCKMATRIX_H 1080