![]() |
Eigen-unsupported
3.3.3
|
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