BlockSparseMatrix.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2013 Desire Nuentsa <desire.nuentsa_wakam@inria.fr>
00005 // Copyright (C) 2013 Gael Guennebaud <gael.guennebaud@inria.fr>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_SPARSEBLOCKMATRIX_H
00012 #define EIGEN_SPARSEBLOCKMATRIX_H
00013 
00014 namespace Eigen { 
00054 template<typename _Scalar, int _BlockAtCompileTime=Dynamic, int _Options=ColMajor, typename _StorageIndex=int> class BlockSparseMatrix;
00055 
00056 template<typename BlockSparseMatrixT> class BlockSparseMatrixView;
00057 
00058 namespace internal {
00059 template<typename _Scalar, int _BlockAtCompileTime, int _Options, typename _Index>
00060 struct traits<BlockSparseMatrix<_Scalar,_BlockAtCompileTime,_Options, _Index> >
00061 {
00062   typedef _Scalar Scalar;
00063   typedef _Index Index;
00064   typedef Sparse StorageKind; // FIXME Where is it used ??
00065   typedef MatrixXpr XprKind;
00066   enum {
00067     RowsAtCompileTime = Dynamic,
00068     ColsAtCompileTime = Dynamic,
00069     MaxRowsAtCompileTime = Dynamic,
00070     MaxColsAtCompileTime = Dynamic,
00071     BlockSize = _BlockAtCompileTime,
00072     Flags = _Options | NestByRefBit | LvalueBit,
00073     CoeffReadCost = NumTraits<Scalar>::ReadCost,
00074     SupportedAccessPatterns = InnerRandomAccessPattern
00075   };
00076 };
00077 template<typename BlockSparseMatrixT>
00078 struct traits<BlockSparseMatrixView<BlockSparseMatrixT> >
00079 {
00080   typedef Ref<Matrix<typename BlockSparseMatrixT::Scalar, BlockSparseMatrixT::BlockSize, BlockSparseMatrixT::BlockSize> > Scalar;
00081   typedef Ref<Matrix<typename BlockSparseMatrixT::RealScalar, BlockSparseMatrixT::BlockSize, BlockSparseMatrixT::BlockSize> > RealScalar;
00082 
00083 };
00084 
00085 // Function object to sort a triplet list
00086 template<typename Iterator, bool IsColMajor>
00087 struct TripletComp
00088 {
00089   typedef typename Iterator::value_type Triplet;
00090   bool operator()(const Triplet& a, const Triplet& b)
00091   { if(IsColMajor)
00092       return ((a.col() == b.col() && a.row() < b.row()) || (a.col() < b.col()));
00093     else
00094       return ((a.row() == b.row() && a.col() < b.col()) || (a.row() < b.row()));
00095   }
00096 };
00097 } // end namespace internal
00098 
00099 
00100 /* Proxy to view the block sparse matrix as a regular sparse matrix */
00101 template<typename BlockSparseMatrixT>
00102 class BlockSparseMatrixView : public SparseMatrixBase<BlockSparseMatrixT>
00103 {
00104   public:
00105     typedef Ref<typename BlockSparseMatrixT::BlockScalar> Scalar;
00106     typedef Ref<typename BlockSparseMatrixT::BlockRealScalar> RealScalar;
00107     typedef typename BlockSparseMatrixT::Index Index;
00108     typedef  BlockSparseMatrixT Nested;
00109     enum {
00110       Flags = BlockSparseMatrixT::Options,
00111       Options = BlockSparseMatrixT::Options,
00112       RowsAtCompileTime = BlockSparseMatrixT::RowsAtCompileTime,
00113       ColsAtCompileTime = BlockSparseMatrixT::ColsAtCompileTime,
00114       MaxColsAtCompileTime = BlockSparseMatrixT::MaxColsAtCompileTime,
00115       MaxRowsAtCompileTime = BlockSparseMatrixT::MaxRowsAtCompileTime
00116     };
00117   public:
00118     BlockSparseMatrixView(const BlockSparseMatrixT& spblockmat)
00119      : m_spblockmat(spblockmat)
00120     {}
00121 
00122     Index outerSize() const
00123     {
00124       return (Flags&RowMajorBit) == 1 ? this->rows() : this->cols();
00125     }
00126     Index cols() const
00127     {
00128       return m_spblockmat.blockCols();
00129     }
00130     Index rows() const
00131     {
00132       return m_spblockmat.blockRows();
00133     }
00134     Scalar coeff(Index row, Index col)
00135     {
00136       return m_spblockmat.coeff(row, col);
00137     }
00138     Scalar coeffRef(Index row, Index col)
00139     {
00140       return m_spblockmat.coeffRef(row, col);
00141     }
00142     // Wrapper to iterate over all blocks
00143     class InnerIterator : public BlockSparseMatrixT::BlockInnerIterator
00144     {
00145       public:
00146       InnerIterator(const BlockSparseMatrixView& mat, Index outer)
00147           : BlockSparseMatrixT::BlockInnerIterator(mat.m_spblockmat, outer)
00148       {}
00149 
00150     };
00151 
00152   protected:
00153     const BlockSparseMatrixT& m_spblockmat;
00154 };
00155 
00156 // Proxy to view a regular vector as a block vector
00157 template<typename BlockSparseMatrixT, typename VectorType>
00158 class BlockVectorView
00159 {
00160   public:
00161     enum {
00162       BlockSize = BlockSparseMatrixT::BlockSize,
00163       ColsAtCompileTime = VectorType::ColsAtCompileTime,
00164       RowsAtCompileTime = VectorType::RowsAtCompileTime,
00165       Flags = VectorType::Flags
00166     };
00167     typedef Ref<const Matrix<typename BlockSparseMatrixT::Scalar, (RowsAtCompileTime==1)? 1 : BlockSize, (ColsAtCompileTime==1)? 1 : BlockSize> >Scalar;
00168     typedef typename BlockSparseMatrixT::Index Index;
00169   public:
00170     BlockVectorView(const BlockSparseMatrixT& spblockmat, const VectorType& vec)
00171     : m_spblockmat(spblockmat),m_vec(vec)
00172     { }
00173     inline Index cols() const
00174     {
00175       return m_vec.cols();
00176     }
00177     inline Index size() const
00178     {
00179       return m_spblockmat.blockRows();
00180     }
00181     inline Scalar coeff(Index bi) const
00182     {
00183       Index startRow = m_spblockmat.blockRowsIndex(bi);
00184       Index rowSize = m_spblockmat.blockRowsIndex(bi+1) - startRow;
00185       return m_vec.middleRows(startRow, rowSize);
00186     }
00187     inline Scalar coeff(Index bi, Index j) const
00188     {
00189       Index startRow = m_spblockmat.blockRowsIndex(bi);
00190       Index rowSize = m_spblockmat.blockRowsIndex(bi+1) - startRow;
00191       return m_vec.block(startRow, j, rowSize, 1);
00192     }
00193   protected:
00194     const BlockSparseMatrixT& m_spblockmat;
00195     const VectorType& m_vec;
00196 };
00197 
00198 template<typename VectorType, typename Index> class BlockVectorReturn;
00199 
00200 
00201 // Proxy to view a regular vector as a block vector
00202 template<typename BlockSparseMatrixT, typename VectorType>
00203 class BlockVectorReturn
00204 {
00205   public:
00206     enum {
00207       ColsAtCompileTime = VectorType::ColsAtCompileTime,
00208       RowsAtCompileTime = VectorType::RowsAtCompileTime,
00209       Flags = VectorType::Flags
00210     };
00211     typedef Ref<Matrix<typename VectorType::Scalar, RowsAtCompileTime, ColsAtCompileTime> > Scalar;
00212     typedef typename BlockSparseMatrixT::Index Index;
00213   public:
00214     BlockVectorReturn(const BlockSparseMatrixT& spblockmat, VectorType& vec)
00215     : m_spblockmat(spblockmat),m_vec(vec)
00216     { }
00217     inline Index size() const
00218     {
00219       return m_spblockmat.blockRows();
00220     }
00221     inline Scalar coeffRef(Index bi)
00222     {
00223       Index startRow = m_spblockmat.blockRowsIndex(bi);
00224       Index rowSize = m_spblockmat.blockRowsIndex(bi+1) - startRow;
00225       return m_vec.middleRows(startRow, rowSize);
00226     }
00227     inline Scalar coeffRef(Index bi, Index j)
00228     {
00229       Index startRow = m_spblockmat.blockRowsIndex(bi);
00230       Index rowSize = m_spblockmat.blockRowsIndex(bi+1) - startRow;
00231       return m_vec.block(startRow, j, rowSize, 1);
00232     }
00233 
00234   protected:
00235     const BlockSparseMatrixT& m_spblockmat;
00236     VectorType& m_vec;
00237 };
00238 
00239 // Block version of the sparse dense product
00240 template<typename Lhs, typename Rhs>
00241 class BlockSparseTimeDenseProduct;
00242 
00243 namespace internal {
00244 
00245 template<typename BlockSparseMatrixT, typename VecType>
00246 struct traits<BlockSparseTimeDenseProduct<BlockSparseMatrixT, VecType> >
00247 {
00248   typedef Dense StorageKind;
00249   typedef MatrixXpr XprKind;
00250   typedef typename BlockSparseMatrixT::Scalar Scalar;
00251   typedef typename BlockSparseMatrixT::Index Index;
00252   enum {
00253     RowsAtCompileTime = Dynamic,
00254     ColsAtCompileTime = Dynamic,
00255     MaxRowsAtCompileTime = Dynamic,
00256     MaxColsAtCompileTime = Dynamic,
00257     Flags = 0,
00258     CoeffReadCost = internal::traits<BlockSparseMatrixT>::CoeffReadCost
00259   };
00260 };
00261 } // end namespace internal
00262 
00263 template<typename Lhs, typename Rhs>
00264 class BlockSparseTimeDenseProduct
00265   : public ProductBase<BlockSparseTimeDenseProduct<Lhs,Rhs>, Lhs, Rhs>
00266 {
00267   public:
00268     EIGEN_PRODUCT_PUBLIC_INTERFACE(BlockSparseTimeDenseProduct)
00269 
00270     BlockSparseTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
00271     {}
00272 
00273     template<typename Dest> void scaleAndAddTo(Dest& dest, const typename Rhs::Scalar& alpha) const
00274     {
00275       BlockVectorReturn<Lhs,Dest> tmpDest(m_lhs, dest);
00276       internal::sparse_time_dense_product( BlockSparseMatrixView<Lhs>(m_lhs),  BlockVectorView<Lhs, Rhs>(m_lhs, m_rhs), tmpDest, alpha);
00277     }
00278 
00279   private:
00280     BlockSparseTimeDenseProduct& operator=(const BlockSparseTimeDenseProduct&);
00281 };
00282 
00283 template<typename _Scalar, int _BlockAtCompileTime, int _Options, typename _StorageIndex>
00284 class BlockSparseMatrix : public SparseMatrixBase<BlockSparseMatrix<_Scalar,_BlockAtCompileTime, _Options,_StorageIndex> >
00285 {
00286   public:
00287     typedef _Scalar Scalar;
00288     typedef typename NumTraits<Scalar>::Real RealScalar;
00289     typedef _StorageIndex StorageIndex;
00290     typedef typename internal::ref_selector<BlockSparseMatrix<_Scalar, _BlockAtCompileTime, _Options, _StorageIndex> >::type Nested;
00291 
00292     enum {
00293       Options = _Options,
00294       Flags = Options,
00295       BlockSize=_BlockAtCompileTime,
00296       RowsAtCompileTime = Dynamic,
00297       ColsAtCompileTime = Dynamic,
00298       MaxRowsAtCompileTime = Dynamic,
00299       MaxColsAtCompileTime = Dynamic,
00300       IsVectorAtCompileTime = 0,
00301       IsColMajor = Flags&RowMajorBit ? 0 : 1
00302     };
00303     typedef Matrix<Scalar, _BlockAtCompileTime, _BlockAtCompileTime,IsColMajor ? ColMajor : RowMajor> BlockScalar;
00304     typedef Matrix<RealScalar, _BlockAtCompileTime, _BlockAtCompileTime,IsColMajor ? ColMajor : RowMajor> BlockRealScalar;
00305     typedef typename internal::conditional<_BlockAtCompileTime==Dynamic, Scalar, BlockScalar>::type BlockScalarReturnType;
00306     typedef BlockSparseMatrix<Scalar, BlockSize, IsColMajor ? ColMajor : RowMajor, StorageIndex> PlainObject;
00307   public:
00308     // Default constructor
00309     BlockSparseMatrix()
00310     : m_innerBSize(0),m_outerBSize(0),m_innerOffset(0),m_outerOffset(0),
00311       m_nonzerosblocks(0),m_values(0),m_blockPtr(0),m_indices(0),
00312       m_outerIndex(0),m_blockSize(BlockSize)
00313     { }
00314 
00315 
00320     BlockSparseMatrix(Index brow, Index bcol)
00321       : m_innerBSize(IsColMajor ? brow : bcol),
00322         m_outerBSize(IsColMajor ? bcol : brow),
00323         m_innerOffset(0),m_outerOffset(0),m_nonzerosblocks(0),
00324         m_values(0),m_blockPtr(0),m_indices(0),
00325         m_outerIndex(0),m_blockSize(BlockSize)
00326     { }
00327 
00331     BlockSparseMatrix(const BlockSparseMatrix& other)
00332       : m_innerBSize(other.m_innerBSize),m_outerBSize(other.m_outerBSize),
00333         m_nonzerosblocks(other.m_nonzerosblocks),m_nonzeros(other.m_nonzeros),
00334         m_blockPtr(0),m_blockSize(other.m_blockSize)
00335     {
00336       // should we allow copying between variable-size blocks and fixed-size blocks ??
00337       eigen_assert(m_blockSize == BlockSize && " CAN NOT COPY BETWEEN FIXED-SIZE AND VARIABLE-SIZE BLOCKS");
00338 
00339       std::copy(other.m_innerOffset, other.m_innerOffset+m_innerBSize+1, m_innerOffset);
00340       std::copy(other.m_outerOffset, other.m_outerOffset+m_outerBSize+1, m_outerOffset);
00341       std::copy(other.m_values, other.m_values+m_nonzeros, m_values);
00342 
00343       if(m_blockSize != Dynamic)
00344         std::copy(other.m_blockPtr, other.m_blockPtr+m_nonzerosblocks, m_blockPtr);
00345 
00346       std::copy(other.m_indices, other.m_indices+m_nonzerosblocks, m_indices);
00347       std::copy(other.m_outerIndex, other.m_outerIndex+m_outerBSize, m_outerIndex);
00348     }
00349 
00350     friend void swap(BlockSparseMatrix& first, BlockSparseMatrix& second)
00351     {
00352       std::swap(first.m_innerBSize, second.m_innerBSize);
00353       std::swap(first.m_outerBSize, second.m_outerBSize);
00354       std::swap(first.m_innerOffset, second.m_innerOffset);
00355       std::swap(first.m_outerOffset, second.m_outerOffset);
00356       std::swap(first.m_nonzerosblocks, second.m_nonzerosblocks);
00357       std::swap(first.m_nonzeros, second.m_nonzeros);
00358       std::swap(first.m_values, second.m_values);
00359       std::swap(first.m_blockPtr, second.m_blockPtr);
00360       std::swap(first.m_indices, second.m_indices);
00361       std::swap(first.m_outerIndex, second.m_outerIndex);
00362       std::swap(first.m_BlockSize, second.m_blockSize);
00363     }
00364 
00365     BlockSparseMatrix& operator=(BlockSparseMatrix other)
00366     {
00367       //Copy-and-swap paradigm ... avoid leaked data if thrown
00368       swap(*this, other);
00369       return *this;
00370     }
00371 
00372     // Destructor
00373     ~BlockSparseMatrix()
00374     {
00375       delete[] m_outerIndex;
00376       delete[] m_innerOffset;
00377       delete[] m_outerOffset;
00378       delete[] m_indices;
00379       delete[] m_blockPtr;
00380       delete[] m_values;
00381     }
00382 
00383 
00388     template<typename MatrixType>
00389     inline BlockSparseMatrix(const MatrixType& spmat) : m_blockSize(BlockSize)
00390     {
00391       EIGEN_STATIC_ASSERT((m_blockSize != Dynamic), THIS_METHOD_IS_ONLY_FOR_FIXED_SIZE);
00392 
00393       *this = spmat;
00394     }
00395 
00404     template<typename MatrixType>
00405     inline BlockSparseMatrix& operator=(const MatrixType& spmat)
00406     {
00407       eigen_assert((m_innerBSize != 0 && m_outerBSize != 0)
00408                    && "Trying to assign to a zero-size matrix, call resize() first");
00409       eigen_assert(((MatrixType::Options&RowMajorBit) != IsColMajor) && "Wrong storage order");
00410       typedef SparseMatrix<bool,MatrixType::Options,typename MatrixType::Index> MatrixPatternType;
00411       MatrixPatternType  blockPattern(blockRows(), blockCols());
00412       m_nonzeros = 0;
00413 
00414       // First, compute the number of nonzero blocks and their locations
00415       for(StorageIndex bj = 0; bj < m_outerBSize; ++bj)
00416       {
00417         // Browse each outer block and compute the structure
00418         std::vector<bool> nzblocksFlag(m_innerBSize,false);  // Record the existing blocks
00419         blockPattern.startVec(bj);
00420         for(StorageIndex j = blockOuterIndex(bj); j < blockOuterIndex(bj+1); ++j)
00421         {
00422           typename MatrixType::InnerIterator it_spmat(spmat, j);
00423           for(; it_spmat; ++it_spmat)
00424           {
00425             StorageIndex bi = innerToBlock(it_spmat.index()); // Index of the current nonzero block
00426             if(!nzblocksFlag[bi])
00427             {
00428               // Save the index of this nonzero block
00429               nzblocksFlag[bi] = true;
00430               blockPattern.insertBackByOuterInnerUnordered(bj, bi) = true;
00431               // Compute the total number of nonzeros (including explicit zeros in blocks)
00432               m_nonzeros += blockOuterSize(bj) * blockInnerSize(bi);
00433             }
00434           }
00435         } // end current outer block
00436       }
00437       blockPattern.finalize();
00438 
00439       // Allocate the internal arrays
00440       setBlockStructure(blockPattern);
00441 
00442       for(StorageIndex nz = 0; nz < m_nonzeros; ++nz) m_values[nz] = Scalar(0);
00443       for(StorageIndex bj = 0; bj < m_outerBSize; ++bj)
00444       {
00445         // Now copy the values
00446         for(StorageIndex j = blockOuterIndex(bj); j < blockOuterIndex(bj+1); ++j)
00447         {
00448           // Browse the outer block column by column (for column-major matrices)
00449           typename MatrixType::InnerIterator it_spmat(spmat, j);
00450           for(; it_spmat; ++it_spmat)
00451           {
00452             StorageIndex idx = 0; // Position of this block in the column block
00453             StorageIndex bi = innerToBlock(it_spmat.index()); // Index of the current nonzero block
00454             // Go to the inner block where this element belongs to
00455             while(bi > m_indices[m_outerIndex[bj]+idx]) ++idx; // Not expensive for ordered blocks
00456             StorageIndex idxVal;// Get the right position in the array of values for this element
00457             if(m_blockSize == Dynamic)
00458             {
00459               // Offset from all blocks before ...
00460               idxVal =  m_blockPtr[m_outerIndex[bj]+idx];
00461               // ... and offset inside the block
00462               idxVal += (j - blockOuterIndex(bj)) * blockOuterSize(bj) + it_spmat.index() - m_innerOffset[bi];
00463             }
00464             else
00465             {
00466               // All blocks before
00467               idxVal = (m_outerIndex[bj] + idx) * m_blockSize * m_blockSize;
00468               // inside the block
00469               idxVal += (j - blockOuterIndex(bj)) * m_blockSize + (it_spmat.index()%m_blockSize);
00470             }
00471             // Insert the value
00472             m_values[idxVal] = it_spmat.value();
00473           } // end of this column
00474         } // end of this block
00475       } // end of this outer block
00476 
00477       return *this;
00478     }
00479 
00497     template<typename MatrixType>
00498     void setBlockStructure(const MatrixType& blockPattern)
00499     {
00500       resize(blockPattern.rows(), blockPattern.cols());
00501       reserve(blockPattern.nonZeros());
00502 
00503       // Browse the block pattern and set up the various pointers
00504       m_outerIndex[0] = 0;
00505       if(m_blockSize == Dynamic) m_blockPtr[0] = 0;
00506       for(StorageIndex nz = 0; nz < m_nonzeros; ++nz) m_values[nz] = Scalar(0);
00507       for(StorageIndex bj = 0; bj < m_outerBSize; ++bj)
00508       {
00509         //Browse each outer block
00510 
00511         //First, copy and save the indices of nonzero blocks
00512         //FIXME : find a way to avoid this ...
00513         std::vector<int> nzBlockIdx;
00514         typename MatrixType::InnerIterator it(blockPattern, bj);
00515         for(; it; ++it)
00516         {
00517           nzBlockIdx.push_back(it.index());
00518         }
00519         std::sort(nzBlockIdx.begin(), nzBlockIdx.end());
00520 
00521         // Now, fill block indices and (eventually) pointers to blocks
00522         for(StorageIndex idx = 0; idx < nzBlockIdx.size(); ++idx)
00523         {
00524           StorageIndex offset = m_outerIndex[bj]+idx; // offset in m_indices
00525           m_indices[offset] = nzBlockIdx[idx];
00526           if(m_blockSize == Dynamic)
00527             m_blockPtr[offset] = m_blockPtr[offset-1] + blockInnerSize(nzBlockIdx[idx]) * blockOuterSize(bj);
00528           // There is no blockPtr for fixed-size blocks... not needed !???
00529         }
00530         // Save the pointer to the next outer block
00531         m_outerIndex[bj+1] = m_outerIndex[bj] + nzBlockIdx.size();
00532       }
00533     }
00534 
00538     inline void resize(Index brow, Index bcol)
00539     {
00540       m_innerBSize = IsColMajor ? brow : bcol;
00541       m_outerBSize = IsColMajor ? bcol : brow;
00542     }
00543 
00549     inline void setBlockSize(Index blockSize)
00550     {
00551       m_blockSize = blockSize;
00552     }
00553 
00563     inline void setBlockLayout(const VectorXi& rowBlocks, const VectorXi& colBlocks)
00564     {
00565       const VectorXi& innerBlocks = IsColMajor ? rowBlocks : colBlocks;
00566       const VectorXi& outerBlocks = IsColMajor ? colBlocks : rowBlocks;
00567       eigen_assert(m_innerBSize == innerBlocks.size() && "CHECK THE NUMBER OF ROW OR COLUMN BLOCKS");
00568       eigen_assert(m_outerBSize == outerBlocks.size() && "CHECK THE NUMBER OF ROW OR COLUMN BLOCKS");
00569       m_outerBSize = outerBlocks.size();
00570       //  starting index of blocks... cumulative sums
00571       m_innerOffset = new StorageIndex[m_innerBSize+1];
00572       m_outerOffset = new StorageIndex[m_outerBSize+1];
00573       m_innerOffset[0] = 0;
00574       m_outerOffset[0] = 0;
00575       std::partial_sum(&innerBlocks[0], &innerBlocks[m_innerBSize-1]+1, &m_innerOffset[1]);
00576       std::partial_sum(&outerBlocks[0], &outerBlocks[m_outerBSize-1]+1, &m_outerOffset[1]);
00577 
00578       // Compute the total number of nonzeros
00579       m_nonzeros = 0;
00580       for(StorageIndex bj = 0; bj < m_outerBSize; ++bj)
00581         for(StorageIndex bi = 0; bi < m_innerBSize; ++bi)
00582           m_nonzeros += outerBlocks[bj] * innerBlocks[bi];
00583 
00584     }
00585 
00596     inline void reserve(const Index nonzerosblocks)
00597     {
00598       eigen_assert((m_innerBSize != 0 && m_outerBSize != 0) &&
00599           "TRYING TO RESERVE ZERO-SIZE MATRICES, CALL resize() first");
00600 
00601       //FIXME Should free if already allocated
00602       m_outerIndex = new StorageIndex[m_outerBSize+1];
00603 
00604       m_nonzerosblocks = nonzerosblocks;
00605       if(m_blockSize != Dynamic)
00606       {
00607         m_nonzeros = nonzerosblocks * (m_blockSize * m_blockSize);
00608         m_blockPtr = 0;
00609       }
00610       else
00611       {
00612         // m_nonzeros  is already computed in setBlockLayout()
00613         m_blockPtr = new StorageIndex[m_nonzerosblocks+1];
00614       }
00615       m_indices = new StorageIndex[m_nonzerosblocks+1];
00616       m_values = new Scalar[m_nonzeros];
00617     }
00618 
00619 
00630     template<typename InputIterator>
00631     void setFromTriplets(const InputIterator& begin, const InputIterator& end)
00632     {
00633       eigen_assert((m_innerBSize!=0 && m_outerBSize !=0) && "ZERO BLOCKS, PLEASE CALL resize() before");
00634 
00635       /* First, sort the triplet list
00636         * FIXME This can be unnecessarily expensive since only the inner indices have to be sorted
00637         * The best approach is like in SparseMatrix::setFromTriplets()
00638         */
00639       internal::TripletComp<InputIterator, IsColMajor> tripletcomp;
00640       std::sort(begin, end, tripletcomp);
00641 
00642       /* Count the number of rows and column blocks,
00643        * and the number of nonzero blocks per outer dimension
00644        */
00645       VectorXi rowBlocks(m_innerBSize); // Size of each block row
00646       VectorXi colBlocks(m_outerBSize); // Size of each block column
00647       rowBlocks.setZero(); colBlocks.setZero();
00648       VectorXi nzblock_outer(m_outerBSize); // Number of nz blocks per outer vector
00649       VectorXi nz_outer(m_outerBSize); // Number of nz per outer vector...for variable-size blocks
00650       nzblock_outer.setZero();
00651       nz_outer.setZero();
00652       for(InputIterator it(begin); it !=end; ++it)
00653       {
00654         eigen_assert(it->row() >= 0 && it->row() < this->blockRows() && it->col() >= 0 && it->col() < this->blockCols());
00655         eigen_assert((it->value().rows() == it->value().cols() && (it->value().rows() == m_blockSize))
00656                      || (m_blockSize == Dynamic));
00657 
00658         if(m_blockSize == Dynamic)
00659         {
00660           eigen_assert((rowBlocks[it->row()] == 0 || rowBlocks[it->row()] == it->value().rows()) &&
00661               "NON CORRESPONDING SIZES FOR ROW BLOCKS");
00662           eigen_assert((colBlocks[it->col()] == 0 || colBlocks[it->col()] == it->value().cols()) &&
00663               "NON CORRESPONDING SIZES FOR COLUMN BLOCKS");
00664           rowBlocks[it->row()] =it->value().rows();
00665           colBlocks[it->col()] = it->value().cols();
00666         }
00667         nz_outer(IsColMajor ? it->col() : it->row()) += it->value().rows() * it->value().cols();
00668         nzblock_outer(IsColMajor ? it->col() : it->row())++;
00669       }
00670       // Allocate member arrays
00671       if(m_blockSize == Dynamic) setBlockLayout(rowBlocks, colBlocks);
00672       StorageIndex nzblocks = nzblock_outer.sum();
00673       reserve(nzblocks);
00674 
00675        // Temporary markers
00676       VectorXi block_id(m_outerBSize); // To be used as a block marker during insertion
00677 
00678       // Setup outer index pointers and markers
00679       m_outerIndex[0] = 0;
00680       if (m_blockSize == Dynamic)  m_blockPtr[0] =  0;
00681       for(StorageIndex bj = 0; bj < m_outerBSize; ++bj)
00682       {
00683         m_outerIndex[bj+1] = m_outerIndex[bj] + nzblock_outer(bj);
00684         block_id(bj) = m_outerIndex[bj];
00685         if(m_blockSize==Dynamic)
00686         {
00687           m_blockPtr[m_outerIndex[bj+1]] = m_blockPtr[m_outerIndex[bj]] + nz_outer(bj);
00688         }
00689       }
00690 
00691       // Fill the matrix
00692       for(InputIterator it(begin); it!=end; ++it)
00693       {
00694         StorageIndex outer = IsColMajor ? it->col() : it->row();
00695         StorageIndex inner = IsColMajor ? it->row() : it->col();
00696         m_indices[block_id(outer)] = inner;
00697         StorageIndex block_size = it->value().rows()*it->value().cols();
00698         StorageIndex nz_marker = blockPtr(block_id[outer]);
00699         memcpy(&(m_values[nz_marker]), it->value().data(), block_size * sizeof(Scalar));
00700         if(m_blockSize == Dynamic)
00701         {
00702           m_blockPtr[block_id(outer)+1] = m_blockPtr[block_id(outer)] + block_size;
00703         }
00704         block_id(outer)++;
00705       }
00706 
00707       // An alternative when the outer indices are sorted...no need to use an array of markers
00708 //      for(Index bcol = 0; bcol < m_outerBSize; ++bcol)
00709 //      {
00710 //      Index id = 0, id_nz = 0, id_nzblock = 0;
00711 //      for(InputIterator it(begin); it!=end; ++it)
00712 //      {
00713 //        while (id<bcol) // one pass should do the job unless there are empty columns
00714 //        {
00715 //          id++;
00716 //          m_outerIndex[id+1]=m_outerIndex[id];
00717 //        }
00718 //        m_outerIndex[id+1] += 1;
00719 //        m_indices[id_nzblock]=brow;
00720 //        Index block_size = it->value().rows()*it->value().cols();
00721 //        m_blockPtr[id_nzblock+1] = m_blockPtr[id_nzblock] + block_size;
00722 //        id_nzblock++;
00723 //        memcpy(&(m_values[id_nz]),it->value().data(), block_size*sizeof(Scalar));
00724 //        id_nz += block_size;
00725 //      }
00726 //      while(id < m_outerBSize-1) // Empty columns at the end
00727 //      {
00728 //        id++;
00729 //        m_outerIndex[id+1]=m_outerIndex[id];
00730 //      }
00731 //      }
00732     }
00733 
00734 
00738     inline Index rows() const
00739     {
00740 //      return blockRows();
00741       return (IsColMajor ? innerSize() : outerSize());
00742     }
00743 
00747     inline Index cols() const
00748     {
00749 //      return blockCols();
00750       return (IsColMajor ? outerSize() : innerSize());
00751     }
00752 
00753     inline Index innerSize() const
00754     {
00755       if(m_blockSize == Dynamic) return m_innerOffset[m_innerBSize];
00756       else return  (m_innerBSize * m_blockSize) ;
00757     }
00758 
00759     inline Index outerSize() const
00760     {
00761       if(m_blockSize == Dynamic) return m_outerOffset[m_outerBSize];
00762       else return  (m_outerBSize * m_blockSize) ;
00763     }
00765     inline Index blockRows() const
00766     {
00767       return (IsColMajor ? m_innerBSize : m_outerBSize);
00768     }
00770     inline Index blockCols() const
00771     {
00772       return (IsColMajor ? m_outerBSize : m_innerBSize);
00773     }
00774 
00775     inline Index outerBlocks() const { return m_outerBSize; }
00776     inline Index innerBlocks() const { return m_innerBSize; }
00777 
00779     inline Index outerToBlock(Index outer) const
00780     {
00781       eigen_assert(outer < outerSize() && "OUTER INDEX OUT OF BOUNDS");
00782 
00783       if(m_blockSize != Dynamic)
00784         return (outer / m_blockSize); // Integer division
00785 
00786       StorageIndex b_outer = 0;
00787       while(m_outerOffset[b_outer] <= outer) ++b_outer;
00788       return b_outer - 1;
00789     }
00791     inline Index innerToBlock(Index inner) const
00792     {
00793       eigen_assert(inner < innerSize() && "OUTER INDEX OUT OF BOUNDS");
00794 
00795       if(m_blockSize != Dynamic)
00796         return (inner / m_blockSize); // Integer division
00797 
00798       StorageIndex b_inner = 0;
00799       while(m_innerOffset[b_inner] <= inner) ++b_inner;
00800       return b_inner - 1;
00801     }
00802 
00806     Ref<BlockScalar> coeffRef(Index brow, Index bcol)
00807     {
00808       eigen_assert(brow < blockRows() && "BLOCK ROW INDEX OUT OF BOUNDS");
00809       eigen_assert(bcol < blockCols() && "BLOCK nzblocksFlagCOLUMN OUT OF BOUNDS");
00810 
00811       StorageIndex rsize = IsColMajor ? blockInnerSize(brow): blockOuterSize(bcol);
00812       StorageIndex csize = IsColMajor ? blockOuterSize(bcol) : blockInnerSize(brow);
00813       StorageIndex inner = IsColMajor ? brow : bcol;
00814       StorageIndex outer = IsColMajor ? bcol : brow;
00815       StorageIndex offset = m_outerIndex[outer];
00816       while(offset < m_outerIndex[outer+1] && m_indices[offset] != inner)
00817         offset++;
00818       if(m_indices[offset] == inner)
00819       {
00820         return Map<BlockScalar>(&(m_values[blockPtr(offset)]), rsize, csize);
00821       }
00822       else
00823       {
00824         //FIXME the block does not exist, Insert it !!!!!!!!!
00825         eigen_assert("DYNAMIC INSERTION IS NOT YET SUPPORTED");
00826       }
00827     }
00828 
00832     Map<const BlockScalar> coeff(Index brow, Index bcol) const
00833     {
00834       eigen_assert(brow < blockRows() && "BLOCK ROW INDEX OUT OF BOUNDS");
00835       eigen_assert(bcol < blockCols() && "BLOCK COLUMN OUT OF BOUNDS");
00836 
00837       StorageIndex rsize = IsColMajor ? blockInnerSize(brow): blockOuterSize(bcol);
00838       StorageIndex csize = IsColMajor ? blockOuterSize(bcol) : blockInnerSize(brow);
00839       StorageIndex inner = IsColMajor ? brow : bcol;
00840       StorageIndex outer = IsColMajor ? bcol : brow;
00841       StorageIndex offset = m_outerIndex[outer];
00842       while(offset < m_outerIndex[outer+1] && m_indices[offset] != inner) offset++;
00843       if(m_indices[offset] == inner)
00844       {
00845         return Map<const BlockScalar> (&(m_values[blockPtr(offset)]), rsize, csize);
00846       }
00847       else
00848 //        return BlockScalar::Zero(rsize, csize);
00849         eigen_assert("NOT YET SUPPORTED");
00850     }
00851 
00852     // Block Matrix times vector product
00853     template<typename VecType>
00854     BlockSparseTimeDenseProduct<BlockSparseMatrix, VecType> operator*(const VecType& lhs) const
00855     {
00856       return BlockSparseTimeDenseProduct<BlockSparseMatrix, VecType>(*this, lhs);
00857     }
00858 
00860     inline Index nonZerosBlocks() const { return m_nonzerosblocks; }
00862     inline Index nonZeros() const { return m_nonzeros; }
00863 
00864     inline BlockScalarReturnType *valuePtr() {return static_cast<BlockScalarReturnType *>(m_values);}
00865 //    inline Scalar *valuePtr(){ return m_values; }
00866     inline StorageIndex *innerIndexPtr() {return m_indices; }
00867     inline const StorageIndex *innerIndexPtr() const {return m_indices; }
00868     inline StorageIndex *outerIndexPtr() {return m_outerIndex; }
00869     inline const StorageIndex* outerIndexPtr() const {return m_outerIndex; }
00870 
00872     inline bool isCompressed() const {return true;}
00876     inline Index blockRowsIndex(Index bi) const
00877     {
00878       return IsColMajor ? blockInnerIndex(bi) : blockOuterIndex(bi);
00879     }
00880 
00884     inline Index blockColsIndex(Index bj) const
00885     {
00886       return IsColMajor ? blockOuterIndex(bj) : blockInnerIndex(bj);
00887     }
00888 
00889     inline Index blockOuterIndex(Index bj) const
00890     {
00891       return (m_blockSize == Dynamic) ? m_outerOffset[bj] : (bj * m_blockSize);
00892     }
00893     inline Index blockInnerIndex(Index bi) const
00894     {
00895       return (m_blockSize == Dynamic) ? m_innerOffset[bi] : (bi * m_blockSize);
00896     }
00897 
00898     // Not needed ???
00899     inline Index blockInnerSize(Index bi) const
00900     {
00901       return (m_blockSize == Dynamic) ? (m_innerOffset[bi+1] - m_innerOffset[bi]) : m_blockSize;
00902     }
00903     inline Index blockOuterSize(Index bj) const
00904     {
00905       return (m_blockSize == Dynamic) ? (m_outerOffset[bj+1]- m_outerOffset[bj]) : m_blockSize;
00906     }
00907 
00911     class InnerIterator; // Browse column by column
00912 
00916     class BlockInnerIterator; // Browse block by block
00917 
00918     friend std::ostream & operator << (std::ostream & s, const BlockSparseMatrix& m)
00919     {
00920       for (StorageIndex j = 0; j < m.outerBlocks(); ++j)
00921       {
00922         BlockInnerIterator itb(m, j);
00923         for(; itb; ++itb)
00924         {
00925           s << "("<<itb.row() << ", " << itb.col() << ")\n";
00926           s << itb.value() <<"\n";
00927         }
00928       }
00929       s << std::endl;
00930       return s;
00931     }
00932 
00936     Index blockPtr(Index id) const
00937     {
00938       if(m_blockSize == Dynamic) return m_blockPtr[id];
00939       else return id * m_blockSize * m_blockSize;
00940       //return blockDynIdx(id, typename internal::conditional<(BlockSize==Dynamic), internal::true_type, internal::false_type>::type());
00941     }
00942 
00943 
00944   protected:
00945 //    inline Index blockDynIdx(Index id, internal::true_type) const
00946 //    {
00947 //      return m_blockPtr[id];
00948 //    }
00949 //    inline Index blockDynIdx(Index id, internal::false_type) const
00950 //    {
00951 //      return id * BlockSize * BlockSize;
00952 //    }
00953 
00954     // To be implemented
00955     // Insert a block at a particular location... need to make a room for that
00956     Map<BlockScalar> insert(Index brow, Index bcol);
00957 
00958     Index m_innerBSize; // Number of block rows
00959     Index m_outerBSize; // Number of block columns
00960     StorageIndex *m_innerOffset; // Starting index of each inner block (size m_innerBSize+1)
00961     StorageIndex *m_outerOffset; // Starting index of each outer block (size m_outerBSize+1)
00962     Index m_nonzerosblocks; // Total nonzeros blocks (lower than  m_innerBSize x m_outerBSize)
00963     Index m_nonzeros; // Total nonzeros elements
00964     Scalar *m_values; //Values stored block column after block column (size m_nonzeros)
00965     StorageIndex *m_blockPtr; // Pointer to the beginning of each block in m_values, size m_nonzeroblocks ... null for fixed-size blocks
00966     StorageIndex *m_indices; //Inner block indices, size m_nonzerosblocks ... OK
00967     StorageIndex *m_outerIndex; // Starting pointer of each block column in m_indices (size m_outerBSize)... OK
00968     Index m_blockSize; // Size of a block for fixed-size blocks, otherwise -1
00969 };
00970 
00971 template<typename _Scalar, int _BlockAtCompileTime, int _Options, typename _StorageIndex>
00972 class BlockSparseMatrix<_Scalar, _BlockAtCompileTime, _Options, _StorageIndex>::BlockInnerIterator
00973 {
00974   public:
00975 
00976     enum{
00977       Flags = _Options
00978     };
00979 
00980     BlockInnerIterator(const BlockSparseMatrix& mat, const Index outer)
00981     : m_mat(mat),m_outer(outer),
00982       m_id(mat.m_outerIndex[outer]),
00983       m_end(mat.m_outerIndex[outer+1])
00984     {
00985     }
00986 
00987     inline BlockInnerIterator& operator++() {m_id++; return *this; }
00988 
00989     inline const Map<const BlockScalar> value() const
00990     {
00991       return Map<const BlockScalar>(&(m_mat.m_values[m_mat.blockPtr(m_id)]),
00992           rows(),cols());
00993     }
00994     inline Map<BlockScalar> valueRef()
00995     {
00996       return Map<BlockScalar>(&(m_mat.m_values[m_mat.blockPtr(m_id)]),
00997           rows(),cols());
00998     }
00999     // Block inner index
01000     inline Index index() const {return m_mat.m_indices[m_id]; }
01001     inline Index outer() const { return m_outer; }
01002     // block row index
01003     inline Index row() const  {return index(); }
01004     // block column index
01005     inline Index col() const {return outer(); }
01006     // FIXME Number of rows in the current block
01007     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; }
01008     // Number of columns in the current block ...
01009     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;}
01010     inline operator bool() const { return (m_id < m_end); }
01011 
01012   protected:
01013     const BlockSparseMatrix<_Scalar, _BlockAtCompileTime, _Options, StorageIndex>& m_mat;
01014     const Index m_outer;
01015     Index m_id;
01016     Index m_end;
01017 };
01018 
01019 template<typename _Scalar, int _BlockAtCompileTime, int _Options, typename _StorageIndex>
01020 class BlockSparseMatrix<_Scalar, _BlockAtCompileTime, _Options, _StorageIndex>::InnerIterator
01021 {
01022   public:
01023     InnerIterator(const BlockSparseMatrix& mat, Index outer)
01024     : m_mat(mat),m_outerB(mat.outerToBlock(outer)),m_outer(outer),
01025       itb(mat, mat.outerToBlock(outer)),
01026       m_offset(outer - mat.blockOuterIndex(m_outerB))
01027      {
01028         if (itb)
01029         {
01030           m_id = m_mat.blockInnerIndex(itb.index());
01031           m_start = m_id;
01032           m_end = m_mat.blockInnerIndex(itb.index()+1);
01033         }
01034      }
01035     inline InnerIterator& operator++()
01036     {
01037       m_id++;
01038       if (m_id >= m_end)
01039       {
01040         ++itb;
01041         if (itb)
01042         {
01043           m_id = m_mat.blockInnerIndex(itb.index());
01044           m_start = m_id;
01045           m_end = m_mat.blockInnerIndex(itb.index()+1);
01046         }
01047       }
01048       return *this;
01049     }
01050     inline const Scalar& value() const
01051     {
01052       return itb.value().coeff(m_id - m_start, m_offset);
01053     }
01054     inline Scalar& valueRef()
01055     {
01056       return itb.valueRef().coeff(m_id - m_start, m_offset);
01057     }
01058     inline Index index() const { return m_id; }
01059     inline Index outer() const {return m_outer; }
01060     inline Index col() const {return outer(); }
01061     inline Index row() const { return index();}
01062     inline operator bool() const
01063     {
01064       return itb;
01065     }
01066   protected:
01067     const BlockSparseMatrix& m_mat;
01068     const Index m_outer;
01069     const Index m_outerB;
01070     BlockInnerIterator itb; // Iterator through the blocks
01071     const Index m_offset; // Position of this column in the block
01072     Index m_start; // starting inner index of this block
01073     Index m_id; // current inner index in the block
01074     Index m_end; // starting inner index of the next block
01075 
01076 };
01077 } // end namespace Eigen
01078 
01079 #endif // EIGEN_SPARSEBLOCKMATRIX_H
 All Classes Functions Variables Typedefs Enumerator