Quelle BlockSparseMatrix.h Sprache: C

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2013 Desire Nuentsa <desire.nuentsa_wakam@inria.fr>
// Copyright (C) 2013 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_SPARSEBLOCKMATRIX_H
#define EIGEN_SPARSEBLOCKMATRIX_H

namespace Eigen {
/** \ingroup SparseCore_Module
  *
  * \class BlockSparseMatrix
  *
  * \brief A versatile sparse matrix representation where each element is a block
  *
  * This class provides routines to manipulate block sparse matrices stored in a
  * BSR-like representation. There are two main types :
  *
  * 1. All blocks have the same number of rows and columns, called block size
  * in the following. In this case, if this block size is known at compile time,
  * it can be given as a template parameter like
  * \code
  * BlockSparseMatrix<Scalar, 3, ColMajor> bmat(b_rows, b_cols);
  * \endcode
  * Here, bmat is a b_rows x b_cols block sparse matrix
  * where each coefficient is a 3x3 dense matrix.
  * If the block size is fixed but will be given at runtime,
  * \code
  * BlockSparseMatrix<Scalar, Dynamic, ColMajor> bmat(b_rows, b_cols);
  * bmat.setBlockSize(block_size);
  * \endcode
  *
  * 2. The second case is for variable-block sparse matrices.
  * Here each block has its own dimensions. The only restriction is that all the blocks
  * in a row (resp. a column) should have the same number of rows (resp. of columns).
  * It is thus required in this case to describe the layout of the matrix by calling
  * setBlockLayout(rowBlocks, colBlocks).
  *
  * In any of the previous case, the matrix can be filled by calling setFromTriplets().
  * A regular sparse matrix can be converted to a block sparse matrix and vice versa.
  * It is obviously required to describe the block layout beforehand by calling either
  * setBlockSize() for fixed-size blocks or setBlockLayout for variable-size blocks.
  *
  * \tparam _Scalar The Scalar type
  * \tparam _BlockAtCompileTime The block layout option. It takes the following values
  * Dynamic : block size known at runtime
  * a numeric number : fixed-size block known at compile time
  */
template<typename _Scalar, int _BlockAtCompileTime=Dynamic, int _Options=ColMajor, typename _StorageIndex=int> class BlockSparseMatrix;

template<typename BlockSparseMatrixT> class BlockSparseMatrixView;

namespace internal {
template<typename _Scalar, int _BlockAtCompileTime, int _Options, typename _Index>
struct traits<BlockSparseMatrix<_Scalar,_BlockAtCompileTime,_Options, _Index> >
{
  typedef _Scalar Scalar;
  typedef _Index Index;
  typedef Sparse StorageKind; // FIXME Where is it used ??
  typedef MatrixXpr XprKind;
  enum {
    RowsAtCompileTime = Dynamic,
    ColsAtCompileTime = Dynamic,
    MaxRowsAtCompileTime = Dynamic,
    MaxColsAtCompileTime = Dynamic,
    BlockSize = _BlockAtCompileTime,
    Flags = _Options | NestByRefBit | LvalueBit,
    CoeffReadCost = NumTraits<Scalar>::ReadCost,
    SupportedAccessPatterns = InnerRandomAccessPattern
  };
};
template<typename BlockSparseMatrixT>
struct traits<BlockSparseMatrixView<BlockSparseMatrixT> >
{
  typedef Ref<Matrix<typename BlockSparseMatrixT::Scalar, BlockSparseMatrixT::BlockSize, BlockSparseMatrixT::BlockSize> > Scalar;
  typedef Ref<Matrix<typename BlockSparseMatrixT::RealScalar, BlockSparseMatrixT::BlockSize, BlockSparseMatrixT::BlockSize> > RealScalar;

};

// Function object to sort a triplet list
template<typename Iterator, bool IsColMajor>
struct TripletComp
{
  typedef typename Iterator::value_type Triplet;
  bool operator()(const Triplet& a, const Triplet& b)
  { if(IsColMajor)
      return ((a.col() == b.col() && a.row() < b.row()) || (a.col() < b.col()));
    else
      return ((a.row() == b.row() && a.col() < b.col()) || (a.row() < b.row()));
  }
};
} // end namespace internal

/* Proxy to view the block sparse matrix as a regular sparse matrix */
template<typename BlockSparseMatrixT>
class BlockSparseMatrixView : public SparseMatrixBase<BlockSparseMatrixT>
{
  public:
    typedef Ref<typename BlockSparseMatrixT::BlockScalar> Scalar;
    typedef Ref<typename BlockSparseMatrixT::BlockRealScalar> RealScalar;
    typedef typename BlockSparseMatrixT::Index Index;
    typedef  BlockSparseMatrixT Nested;
    enum {
      Flags = BlockSparseMatrixT::Options,
      Options = BlockSparseMatrixT::Options,
      RowsAtCompileTime = BlockSparseMatrixT::RowsAtCompileTime,
      ColsAtCompileTime = BlockSparseMatrixT::ColsAtCompileTime,
      MaxColsAtCompileTime = BlockSparseMatrixT::MaxColsAtCompileTime,
      MaxRowsAtCompileTime = BlockSparseMatrixT::MaxRowsAtCompileTime
    };
  public:
    BlockSparseMatrixView(const BlockSparseMatrixT& spblockmat)
     : m_spblockmat(spblockmat)
    {}

    Index outerSize() const
    {
      return (Flags&RowMajorBit) == 1 ? this->rows() : this->cols();
    }
    Index cols() const
    {
      return m_spblockmat.blockCols();
    }
    Index rows() const
    {
      return m_spblockmat.blockRows();
    }
    Scalar coeff(Index row, Index col)
    {
      return m_spblockmat.coeff(row, col);
    }
    Scalar coeffRef(Index row, Index col)
    {
      return m_spblockmat.coeffRef(row, col);
    }
    // Wrapper to iterate over all blocks
    class InnerIterator : public BlockSparseMatrixT::BlockInnerIterator
    {
      public:
      InnerIterator(const BlockSparseMatrixView& mat, Index outer)
          : BlockSparseMatrixT::BlockInnerIterator(mat.m_spblockmat, outer)
      {}

    };

  protected:
    const BlockSparseMatrixT& m_spblockmat;
};

// Proxy to view a regular vector as a block vector
template<typename BlockSparseMatrixT, typename VectorType>
class BlockVectorView
{
  public:
    enum {
      BlockSize = BlockSparseMatrixT::BlockSize,
      ColsAtCompileTime = VectorType::ColsAtCompileTime,
      RowsAtCompileTime = VectorType::RowsAtCompileTime,
      Flags = VectorType::Flags
    };
    typedef Ref<const Matrix<typename BlockSparseMatrixT::Scalar, (RowsAtCompileTime==1)? 1 : BlockSize, (ColsAtCompileTime==1)? 1 : BlockSize> >Scalar;
    typedef typename BlockSparseMatrixT::Index Index;
  public:
    BlockVectorView(const BlockSparseMatrixT& spblockmat, const VectorType& vec)
    : m_spblockmat(spblockmat),m_vec(vec)
    { }
    inline Index cols() const
    {
      return m_vec.cols();
    }
    inline Index size() const
    {
      return m_spblockmat.blockRows();
    }
    inline Scalar coeff(Index bi) const
    {
      Index startRow = m_spblockmat.blockRowsIndex(bi);
      Index rowSize = m_spblockmat.blockRowsIndex(bi+1) - startRow;
      return m_vec.middleRows(startRow, rowSize);
    }
    inline Scalar coeff(Index bi, Index j) const
    {
      Index startRow = m_spblockmat.blockRowsIndex(bi);
      Index rowSize = m_spblockmat.blockRowsIndex(bi+1) - startRow;
      return m_vec.block(startRow, j, rowSize, 1);
    }
  protected:
    const BlockSparseMatrixT& m_spblockmat;
    const VectorType& m_vec;
};

template<typename VectorType, typename Index> class BlockVectorReturn;

// Proxy to view a regular vector as a block vector
template<typename BlockSparseMatrixT, typename VectorType>
class BlockVectorReturn
{
  public:
    enum {
      ColsAtCompileTime = VectorType::ColsAtCompileTime,
      RowsAtCompileTime = VectorType::RowsAtCompileTime,
      Flags = VectorType::Flags
    };
    typedef Ref<Matrix<typename VectorType::Scalar, RowsAtCompileTime, ColsAtCompileTime> > Scalar;
    typedef typename BlockSparseMatrixT::Index Index;
  public:
    BlockVectorReturn(const BlockSparseMatrixT& spblockmat, VectorType& vec)
    : m_spblockmat(spblockmat),m_vec(vec)
    { }
    inline Index size() const
    {
      return m_spblockmat.blockRows();
    }
    inline Scalar coeffRef(Index bi)
    {
      Index startRow = m_spblockmat.blockRowsIndex(bi);
      Index rowSize = m_spblockmat.blockRowsIndex(bi+1) - startRow;
      return m_vec.middleRows(startRow, rowSize);
    }
    inline Scalar coeffRef(Index bi, Index j)
    {
      Index startRow = m_spblockmat.blockRowsIndex(bi);
      Index rowSize = m_spblockmat.blockRowsIndex(bi+1) - startRow;
      return m_vec.block(startRow, j, rowSize, 1);
    }

  protected:
    const BlockSparseMatrixT& m_spblockmat;
    VectorType& m_vec;
};

// Block version of the sparse dense product
template<typename Lhs, typename Rhs>
class BlockSparseTimeDenseProduct;

namespace internal {

template<typename BlockSparseMatrixT, typename VecType>
struct traits<BlockSparseTimeDenseProduct<BlockSparseMatrixT, VecType> >
{
  typedef Dense StorageKind;
  typedef MatrixXpr XprKind;
  typedef typename BlockSparseMatrixT::Scalar Scalar;
  typedef typename BlockSparseMatrixT::Index Index;
  enum {
    RowsAtCompileTime = Dynamic,
    ColsAtCompileTime = Dynamic,
    MaxRowsAtCompileTime = Dynamic,
    MaxColsAtCompileTime = Dynamic,
    Flags = 0,
    CoeffReadCost = internal::traits<BlockSparseMatrixT>::CoeffReadCost
  };
};
} // end namespace internal

template<typename Lhs, typename Rhs>
class BlockSparseTimeDenseProduct
  : public ProductBase<BlockSparseTimeDenseProduct<Lhs,Rhs>, Lhs, Rhs>
{
  public:
    EIGEN_PRODUCT_PUBLIC_INTERFACE(BlockSparseTimeDenseProduct)

    BlockSparseTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
    {}

    template<typename Dest> void scaleAndAddTo(Dest& dest, const typename Rhs::Scalar& ;alpha) const
    {
      BlockVectorReturn<Lhs,Dest> tmpDest(m_lhs, dest);
      internal::sparse_time_dense_product( BlockSparseMatrixView<Lhs>(m_lhs),  BlockVectorView<Lhs, Rhs>(m_lhs, m_rhs), tmpDest, alpha);
    }

  private:
    BlockSparseTimeDenseProduct& operator=(const BlockSparseTimeDenseProduct&);
};

template<typename _Scalar, int _BlockAtCompileTime, int _Options, typename _StorageIndex>
class BlockSparseMatrix : public SparseMatrixBase<BlockSparseMatrix<_Scalar,_BlockAtCompileTime, _Options,_StorageIndex> >
{
  public:
    typedef _Scalar Scalar;
    typedef typename NumTraits<Scalar>::Real RealScalar;
    typedef _StorageIndex StorageIndex;
    typedef typename internal::ref_selector<BlockSparseMatrix<_Scalar, _BlockAtCompileTime, _Options, _StorageIndex> >::type Nested;

    enum {
      Options = _Options,
      Flags = Options,
      BlockSize=_BlockAtCompileTime,
      RowsAtCompileTime = Dynamic,
      ColsAtCompileTime = Dynamic,
      MaxRowsAtCompileTime = Dynamic,
      MaxColsAtCompileTime = Dynamic,
      IsVectorAtCompileTime = 0,
      IsColMajor = Flags&RowMajorBit ? 0 : 1
    };
    typedef Matrix<Scalar, _BlockAtCompileTime, _BlockAtCompileTime,IsColMajor ? ColMajor : RowMajor> BlockScalar;
    typedef Matrix<RealScalar, _BlockAtCompileTime, _BlockAtCompileTime,IsColMajor ? ColMajor : RowMajor> BlockRealScalar;
    typedef typename internal::conditional<_BlockAtCompileTime==Dynamic, Scalar, BlockScalar>::type BlockScalarReturnType;
    typedef BlockSparseMatrix<Scalar, BlockSize, IsColMajor ? ColMajor : RowMajor, StorageIndex> PlainObject;
  public:
    // Default constructor
    BlockSparseMatrix()
    : m_innerBSize(0),m_outerBSize(0),m_innerOffset(0),m_outerOffset(0),
      m_nonzerosblocks(0),m_values(0),m_blockPtr(0),m_indices(0),
      m_outerIndex(0),m_blockSize(BlockSize)
    { }

    /**
     * \brief Construct and resize
     *
     */
    BlockSparseMatrix(Index brow, Index bcol)
      : m_innerBSize(IsColMajor ? brow : bcol),
        m_outerBSize(IsColMajor ? bcol : brow),
        m_innerOffset(0),m_outerOffset(0),m_nonzerosblocks(0),
        m_values(0),m_blockPtr(0),m_indices(0),
        m_outerIndex(0),m_blockSize(BlockSize)
    { }

    /**
     * \brief Copy-constructor
     */
    BlockSparseMatrix(const BlockSparseMatrix& other)
      : m_innerBSize(other.m_innerBSize),m_outerBSize(other.m_outerBSize),
        m_nonzerosblocks(other.m_nonzerosblocks),m_nonzeros(other.m_nonzeros),
        m_blockPtr(0),m_blockSize(other.m_blockSize)
    {
      // should we allow copying between variable-size blocks and fixed-size blocks ??
      eigen_assert(m_blockSize == BlockSize && " CAN NOT COPY BETWEEN FIXED-SIZE AND VARIABLE-SIZE BLOCKS");

      std::copy(other.m_innerOffset, other.m_innerOffset+m_innerBSize+1, m_innerOffset);
      std::copy(other.m_outerOffset, other.m_outerOffset+m_outerBSize+1, m_outerOffset);
      std::copy(other.m_values, other.m_values+m_nonzeros, m_values);

      if(m_blockSize != Dynamic)
        std::copy(other.m_blockPtr, other.m_blockPtr+m_nonzerosblocks, m_blockPtr);

      std::copy(other.m_indices, other.m_indices+m_nonzerosblocks, m_indices);
      std::copy(other.m_outerIndex, other.m_outerIndex+m_outerBSize, m_outerIndex);
    }

    friend void swap(BlockSparseMatrix& first, BlockSparseMatrix& second)
    {
      std::swap(first.m_innerBSize, second.m_innerBSize);
      std::swap(first.m_outerBSize, second.m_outerBSize);
      std::swap(first.m_innerOffset, second.m_innerOffset);
      std::swap(first.m_outerOffset, second.m_outerOffset);
      std::swap(first.m_nonzerosblocks, second.m_nonzerosblocks);
      std::swap(first.m_nonzeros, second.m_nonzeros);
      std::swap(first.m_values, second.m_values);
      std::swap(first.m_blockPtr, second.m_blockPtr);
      std::swap(first.m_indices, second.m_indices);
      std::swap(first.m_outerIndex, second.m_outerIndex);
      std::swap(first.m_BlockSize, second.m_blockSize);
    }

    BlockSparseMatrix& operator=(BlockSparseMatrix other)
    {
      //Copy-and-swap paradigm ... avoid leaked data if thrown
      swap(*this, other);
      return *this;
    }

    // Destructor
    ~BlockSparseMatrix()
    {
      delete[] m_outerIndex;
      delete[] m_innerOffset;
      delete[] m_outerOffset;
      delete[] m_indices;
      delete[] m_blockPtr;
      delete[] m_values;
    }

    /**
      * \brief Constructor from a sparse matrix
      *
      */
    template<typename MatrixType>
    inline BlockSparseMatrix(const MatrixType& spmat) : m_blockSize(BlockSize)
    {
      EIGEN_STATIC_ASSERT((m_blockSize != Dynamic), THIS_METHOD_IS_ONLY_FOR_FIXED_SIZE);

      *this = spmat;
    }

    /**
      * \brief Assignment from a sparse matrix with the same storage order
      *
      * Convert from a sparse matrix to block sparse matrix.
      * \warning Before calling this function, tt is necessary to call
      * either setBlockLayout() (matrices with variable-size blocks)
      * or setBlockSize() (for fixed-size blocks).
      */
    template<typename MatrixType>
    inline BlockSparseMatrix& operator=(const MatrixType& spmat)
    {
      eigen_assert((m_innerBSize != 0 && m_outerBSize != 0)
                   && "Trying to assign to a zero-size matrix, call resize() first");
      eigen_assert(((MatrixType::Options&RowMajorBit) != IsColMajor) && "Wrong storage order");
      typedef SparseMatrix<bool,MatrixType::Options,typename MatrixType::Index> MatrixPatternType;
      MatrixPatternType  blockPattern(blockRows(), blockCols());
      m_nonzeros = 0;

      // First, compute the number of nonzero blocks and their locations
      for(StorageIndex bj = 0; bj < m_outerBSize; ++bj)
      {
        // Browse each outer block and compute the structure
        std::vector<bool> nzblocksFlag(m_innerBSize,false);  // Record the existing blocks
        blockPattern.startVec(bj);
        for(StorageIndex j = blockOuterIndex(bj); j < blockOuterIndex(bj+1); ++j)
        {
          typename MatrixType::InnerIterator it_spmat(spmat, j);
          for(; it_spmat; ++it_spmat)
          {
            StorageIndex bi = innerToBlock(it_spmat.index()); // Index of the current nonzero block
            if(!nzblocksFlag[bi])
            {
              // Save the index of this nonzero block
              nzblocksFlag[bi] = true;
              blockPattern.insertBackByOuterInnerUnordered(bj, bi) = true;
              // Compute the total number of nonzeros (including explicit zeros in blocks)
              m_nonzeros += blockOuterSize(bj) * blockInnerSize(bi);
            }
          }
        } // end current outer block
      }
      blockPattern.finalize();

      // Allocate the internal arrays
      setBlockStructure(blockPattern);

      for(StorageIndex nz = 0; nz < m_nonzeros; ++nz) m_values[nz] = Scalar(0);
      for(StorageIndex bj = 0; bj < m_outerBSize; ++bj)
      {
        // Now copy the values
        for(StorageIndex j = blockOuterIndex(bj); j < blockOuterIndex(bj+1); ++j)
        {
          // Browse the outer block column by column (for column-major matrices)
          typename MatrixType::InnerIterator it_spmat(spmat, j);
          for(; it_spmat; ++it_spmat)
          {
            StorageIndex idx = 0; // Position of this block in the column block
            StorageIndex bi = innerToBlock(it_spmat.index()); // Index of the current nonzero block
            // Go to the inner block where this element belongs to
            while(bi > m_indices[m_outerIndex[bj]+idx]) ++idx; // Not expensive for ordered blocks
            StorageIndex idxVal;// Get the right position in the array of values for this element
            if(m_blockSize == Dynamic)
            {
              // Offset from all blocks before ...
              idxVal =  m_blockPtr[m_outerIndex[bj]+idx];
              // ... and offset inside the block
              idxVal += (j - blockOuterIndex(bj)) * blockOuterSize(bj) + it_spmat.index() - m_innerOffset[bi];
            }
            else
            {
              // All blocks before
              idxVal = (m_outerIndex[bj] + idx) * m_blockSize * m_blockSize;
              // inside the block
              idxVal += (j - blockOuterIndex(bj)) * m_blockSize + (it_spmat.index()%m_blockSize);
            }
            // Insert the value
            m_values[idxVal] = it_spmat.value();
          } // end of this column
        } // end of this block
      } // end of this outer block

      return *this;
    }

    /**
      * \brief Set the nonzero block pattern of the matrix
      *
      * Given a sparse matrix describing the nonzero block pattern,
      * this function prepares the internal pointers for values.
      * After calling this function, any *nonzero* block (bi, bj) can be set
      * with a simple call to coeffRef(bi,bj).
      *
      *
      * \warning Before calling this function, tt is necessary to call
      * either setBlockLayout() (matrices with variable-size blocks)
      * or setBlockSize() (for fixed-size blocks).
      *
      * \param blockPattern Sparse matrix of boolean elements describing the block structure
      *
      * \sa setBlockLayout() \sa setBlockSize()
      */
    template<typename MatrixType>
    void setBlockStructure(const MatrixType& blockPattern)
    {
      resize(blockPattern.rows(), blockPattern.cols());
      reserve(blockPattern.nonZeros());

      // Browse the block pattern and set up the various pointers
      m_outerIndex[0] = 0;
      if(m_blockSize == Dynamic) m_blockPtr[0] = 0;
      for(StorageIndex nz = 0; nz < m_nonzeros; ++nz) m_values[nz] = Scalar(0);
      for(StorageIndex bj = 0; bj < m_outerBSize; ++bj)
      {
        //Browse each outer block

        //First, copy and save the indices of nonzero blocks
        //FIXME : find a way to avoid this ...
        std::vector<int> nzBlockIdx;
        typename MatrixType::InnerIterator it(blockPattern, bj);
        for(; it; ++it)
        {
          nzBlockIdx.push_back(it.index());
        }
        std::sort(nzBlockIdx.begin(), nzBlockIdx.end());

        // Now, fill block indices and (eventually) pointers to blocks
        for(StorageIndex idx = 0; idx < nzBlockIdx.size(); ++idx)
        {
          StorageIndex offset = m_outerIndex[bj]+idx; // offset in m_indices
          m_indices[offset] = nzBlockIdx[idx];
          if(m_blockSize == Dynamic)
            m_blockPtr[offset] = m_blockPtr[offset-1] + blockInnerSize(nzBlockIdx[idx]) * blockOuterSize(bj);
          // There is no blockPtr for fixed-size blocks... not needed !???
        }
        // Save the pointer to the next outer block
        m_outerIndex[bj+1] = m_outerIndex[bj] + nzBlockIdx.size();
      }
    }

    /**
      * \brief Set the number of rows and columns blocks
      */
    inline void resize(Index brow, Index bcol)
    {
      m_innerBSize = IsColMajor ? brow : bcol;
      m_outerBSize = IsColMajor ? bcol : brow;
    }

    /**
      * \brief set the block size at runtime for fixed-size block layout
      *
      * Call this only for fixed-size blocks
      */
    inline void setBlockSize(Index blockSize)
    {
      m_blockSize = blockSize;
    }

    /**
      * \brief Set the row and column block layouts,
      *
      * This function set the size of each row and column block.
      * So this function should be used only for blocks with variable size.
      * \param rowBlocks : Number of rows per row block
      * \param colBlocks : Number of columns per column block
      * \sa resize(), setBlockSize()
      */
    inline void setBlockLayout(const VectorXi& rowBlocks, const VectorXi& colBlocks)
    {
      const VectorXi& innerBlocks = IsColMajor ? rowBlocks : colBlocks;
      const VectorXi& outerBlocks = IsColMajor ? colBlocks : rowBlocks;
      eigen_assert(m_innerBSize == innerBlocks.size() && "CHECK THE NUMBER OF ROW OR COLUMN BLOCKS");
      eigen_assert(m_outerBSize == outerBlocks.size() && "CHECK THE NUMBER OF ROW OR COLUMN BLOCKS");
      m_outerBSize = outerBlocks.size();
      //  starting index of blocks... cumulative sums
      m_innerOffset = new StorageIndex[m_innerBSize+1];
      m_outerOffset = new StorageIndex[m_outerBSize+1];
      m_innerOffset[0] = 0;
      m_outerOffset[0] = 0;
      std::partial_sum(&innerBlocks[0], &innerBlocks[m_innerBSize-1]+1, &m_innerOffset[1]);
      std::partial_sum(&outerBlocks[0], &outerBlocks[m_outerBSize-1]+1, &m_outerOffset[1]);

      // Compute the total number of nonzeros
      m_nonzeros = 0;
      for(StorageIndex bj = 0; bj < m_outerBSize; ++bj)
        for(StorageIndex bi = 0; bi < m_innerBSize; ++bi)
          m_nonzeros += outerBlocks[bj] * innerBlocks[bi];

    }

    /**
      * \brief Allocate the internal array of pointers to blocks and their inner indices
      *
      * \note For fixed-size blocks, call setBlockSize() to set the block.
      * And For variable-size blocks, call setBlockLayout() before using this function
      *
      * \param nonzerosblocks Number of nonzero blocks. The total number of nonzeros is
      * is computed in setBlockLayout() for variable-size blocks
      * \sa setBlockSize()
      */
    inline void reserve(const Index nonzerosblocks)
    {
      eigen_assert((m_innerBSize != 0 && m_outerBSize != 0) &&
          "TRYING TO RESERVE ZERO-SIZE MATRICES, CALL resize() first");

      //FIXME Should free if already allocated
      m_outerIndex = new StorageIndex[m_outerBSize+1];

      m_nonzerosblocks = nonzerosblocks;
      if(m_blockSize != Dynamic)
      {
        m_nonzeros = nonzerosblocks * (m_blockSize * m_blockSize);
        m_blockPtr = 0;
      }
      else
      {
        // m_nonzeros  is already computed in setBlockLayout()
        m_blockPtr = new StorageIndex[m_nonzerosblocks+1];
      }
      m_indices = new StorageIndex[m_nonzerosblocks+1];
      m_values = new Scalar[m_nonzeros];
    }

    /**
      * \brief Fill values in a matrix  from a triplet list.
      *
      * Each triplet item has a block stored in an Eigen dense matrix.
      * The InputIterator class should provide the functions row(), col() and value()
      *
      * \note For fixed-size blocks, call setBlockSize() before this function.
      *
      * FIXME Do not accept duplicates
      */
    template<typename InputIterator>
    void setFromTriplets(const InputIterator& begin, const InputIterator& end)
    {
      eigen_assert((m_innerBSize!=0 && m_outerBSize !=0) && "ZERO BLOCKS, PLEASE CALL resize() before");

      /* First, sort the triplet list
        * FIXME This can be unnecessarily expensive since only the inner indices have to be sorted
        * The best approach is like in SparseMatrix::setFromTriplets()
        */
      internal::TripletComp<InputIterator, IsColMajor> tripletcomp;
      std::sort(begin, end, tripletcomp);

      /* Count the number of rows and column blocks,
       * and the number of nonzero blocks per outer dimension
       */
      VectorXi rowBlocks(m_innerBSize); // Size of each block row
      VectorXi colBlocks(m_outerBSize); // Size of each block column
      rowBlocks.setZero(); colBlocks.setZero();
      VectorXi nzblock_outer(m_outerBSize); // Number of nz blocks per outer vector
      VectorXi nz_outer(m_outerBSize); // Number of nz per outer vector...for variable-size blocks
      nzblock_outer.setZero();
      nz_outer.setZero();
      for(InputIterator it(begin); it !=end; ++it)
      {
        eigen_assert(it->row() >= 0 && it->row() < this->blockRows() && it->col() >= 0 && it->col() < this->blockCols());
        eigen_assert((it->value().rows() == it->value().cols() && (it->value().rows() == m_blockSize))
                     || (m_blockSize == Dynamic));

        if(m_blockSize == Dynamic)
        {
          eigen_assert((rowBlocks[it->row()] == 0 || rowBlocks[it->row()] == it->value().rows()) &&
              "NON CORRESPONDING SIZES FOR ROW BLOCKS");
          eigen_assert((colBlocks[it->col()] == 0 || colBlocks[it->col()] == it->value().cols()) &&
              "NON CORRESPONDING SIZES FOR COLUMN BLOCKS");
          rowBlocks[it->row()] =it->value().rows();
          colBlocks[it->col()] = it->value().cols();
        }
        nz_outer(IsColMajor ? it->col() : it->row()) += it->value().rows() * it->value().cols();
        nzblock_outer(IsColMajor ? it->col() : it->row())++;
      }
      // Allocate member arrays
      if(m_blockSize == Dynamic) setBlockLayout(rowBlocks, colBlocks);
      StorageIndex nzblocks = nzblock_outer.sum();
      reserve(nzblocks);

       // Temporary markers
      VectorXi block_id(m_outerBSize); // To be used as a block marker during insertion

      // Setup outer index pointers and markers
      m_outerIndex[0] = 0;
      if (m_blockSize == Dynamic)  m_blockPtr[0] =  0;
      for(StorageIndex bj = 0; bj < m_outerBSize; ++bj)
      {
        m_outerIndex[bj+1] = m_outerIndex[bj] + nzblock_outer(bj);
        block_id(bj) = m_outerIndex[bj];
        if(m_blockSize==Dynamic)
        {
          m_blockPtr[m_outerIndex[bj+1]] = m_blockPtr[m_outerIndex[bj]] + nz_outer(bj);
        }
      }

      // Fill the matrix
      for(InputIterator it(begin); it!=end; ++it)
      {
        StorageIndex outer = IsColMajor ? it->col() : it->row();
        StorageIndex inner = IsColMajor ? it->row() : it->col();
        m_indices[block_id(outer)] = inner;
        StorageIndex block_size = it->value().rows()*it->value().cols();
        StorageIndex nz_marker = blockPtr(block_id[outer]);
        memcpy(&(m_values[nz_marker]), it->value().data(), block_size * sizeof(Scalar));
        if(m_blockSize == Dynamic)
        {
          m_blockPtr[block_id(outer)+1] = m_blockPtr[block_id(outer)] + block_size;
        }
        block_id(outer)++;
      }

      // An alternative when the outer indices are sorted...no need to use an array of markers
//      for(Index bcol = 0; bcol < m_outerBSize; ++bcol)
//      {
//      Index id = 0, id_nz = 0, id_nzblock = 0;
//      for(InputIterator it(begin); it!=end; ++it)
//      {
//        while (id<bcol) // one pass should do the job unless there are empty columns
//        {
//          id++;
//          m_outerIndex[id+1]=m_outerIndex[id];
//        }
//        m_outerIndex[id+1] += 1;
//        m_indices[id_nzblock]=brow;
//        Index block_size = it->value().rows()*it->value().cols();
//        m_blockPtr[id_nzblock+1] = m_blockPtr[id_nzblock] + block_size;
//        id_nzblock++;
//        memcpy(&(m_values[id_nz]),it->value().data(), block_size*sizeof(Scalar));
//        id_nz += block_size;
//      }
//      while(id < m_outerBSize-1) // Empty columns at the end
//      {
//        id++;
//        m_outerIndex[id+1]=m_outerIndex[id];
//      }
//      }
    }

    /**
      * \returns the number of rows
      */
    inline Index rows() const
    {
//      return blockRows();
      return (IsColMajor ? innerSize() : outerSize());
    }

    /**
      * \returns the number of cols
      */
    inline Index cols() const
    {
//      return blockCols();
      return (IsColMajor ? outerSize() : innerSize());
    }

    inline Index innerSize() const
    {
      if(m_blockSize == Dynamic) return m_innerOffset[m_innerBSize];
      else return  (m_innerBSize * m_blockSize) ;
    }

    inline Index outerSize() const
    {
      if(m_blockSize == Dynamic) return m_outerOffset[m_outerBSize];
      else return  (m_outerBSize * m_blockSize) ;
    }
    /** \returns the number of rows grouped by blocks */
    inline Index blockRows() const
    {
      return (IsColMajor ? m_innerBSize : m_outerBSize);
    }
    /** \returns the number of columns grouped by blocks */
    inline Index blockCols() const
    {
      return (IsColMajor ? m_outerBSize : m_innerBSize);
    }

    inline Index outerBlocks() const { return m_outerBSize; }
    inline Index innerBlocks() const { return m_innerBSize; }

    /** \returns the block index where outer belongs to */
    inline Index outerToBlock(Index outer) const
    {
      eigen_assert(outer < outerSize() && "OUTER INDEX OUT OF BOUNDS");

      if(m_blockSize != Dynamic)
        return (outer / m_blockSize); // Integer division

      StorageIndex b_outer = 0;
      while(m_outerOffset[b_outer] <= outer) ++b_outer;
      return b_outer - 1;
    }
    /** \returns  the block index where inner belongs to */
    inline Index innerToBlock(Index inner) const
    {
      eigen_assert(inner < innerSize() && "OUTER INDEX OUT OF BOUNDS");

      if(m_blockSize != Dynamic)
        return (inner / m_blockSize); // Integer division

      StorageIndex b_inner = 0;
      while(m_innerOffset[b_inner] <= inner) ++b_inner;
      return b_inner - 1;
    }

    /**
      *\returns a reference to the (i,j) block as an Eigen Dense Matrix
      */
    Ref<BlockScalar> coeffRef(Index brow, Index bcol)
    {
      eigen_assert(brow < blockRows() && "BLOCK ROW INDEX OUT OF BOUNDS");
      eigen_assert(bcol < blockCols() && "BLOCK nzblocksFlagCOLUMN OUT OF BOUNDS");

      StorageIndex rsize = IsColMajor ? blockInnerSize(brow): blockOuterSize(bcol);
      StorageIndex csize = IsColMajor ? blockOuterSize(bcol) : blockInnerSize(brow);
      StorageIndex inner = IsColMajor ? brow : bcol;
      StorageIndex outer = IsColMajor ? bcol : brow;
      StorageIndex offset = m_outerIndex[outer];
      while(offset < m_outerIndex[outer+1] && m_indices[offset] != inner)
        offset++;
      if(m_indices[offset] == inner)
      {
        return Map<BlockScalar>(&(m_values[blockPtr(offset)]), rsize, csize);
      }
      else
      {
        //FIXME the block does not exist, Insert it !!!!!!!!!
        eigen_assert("DYNAMIC INSERTION IS NOT YET SUPPORTED");
      }
    }

    /**
      * \returns the value of the (i,j) block as an Eigen Dense Matrix
      */
    Map<const BlockScalar> coeff(Index brow, Index bcol) const
    {
      eigen_assert(brow < blockRows() && "BLOCK ROW INDEX OUT OF BOUNDS");
      eigen_assert(bcol < blockCols() && "BLOCK COLUMN OUT OF BOUNDS");

      StorageIndex rsize = IsColMajor ? blockInnerSize(brow): blockOuterSize(bcol);
      StorageIndex csize = IsColMajor ? blockOuterSize(bcol) : blockInnerSize(brow);
      StorageIndex inner = IsColMajor ? brow : bcol;
      StorageIndex outer = IsColMajor ? bcol : brow;
      StorageIndex offset = m_outerIndex[outer];
      while(offset < m_outerIndex[outer+1] && m_indices[offset] != inner) offset++;
      if(m_indices[offset] == inner)
      {
        return Map<const BlockScalar> (&(m_values[blockPtr(offset)]), rsize, csize);
      }
      else
//        return BlockScalar::Zero(rsize, csize);
        eigen_assert("NOT YET SUPPORTED");
    }

    // Block Matrix times vector product
    template<typename VecType>
    BlockSparseTimeDenseProduct<BlockSparseMatrix, VecType> operator*(const VecType& lhs) const
    {
      return BlockSparseTimeDenseProduct<BlockSparseMatrix, VecType>(*this, lhs);
    }

    /** \returns the number of nonzero blocks */
    inline Index nonZerosBlocks() const { return m_nonzerosblocks; }
    /** \returns the total number of nonzero elements, including eventual explicit zeros in blocks */
    inline Index nonZeros() const { return m_nonzeros; }

    inline BlockScalarReturnType *valuePtr() {return static_cast<BlockScalarReturnType *>(m_values);}
//    inline Scalar *valuePtr(){ return m_values; }
    inline StorageIndex *innerIndexPtr() {return m_indices; }
    inline const StorageIndex *innerIndexPtr() const {return m_indices; }
    inline StorageIndex *outerIndexPtr() {return m_outerIndex; }
    inline const StorageIndex* outerIndexPtr() const {return m_outerIndex; }

    /** \brief for compatibility purposes with the SparseMatrix class */
    inline bool isCompressed() const {return true;}
    /**
      * \returns the starting index of the bi row block
      */
    inline Index blockRowsIndex(Index bi) const
    {
      return IsColMajor ? blockInnerIndex(bi) : blockOuterIndex(bi);
    }

    /**
      * \returns the starting index of the bj col block
      */
    inline Index blockColsIndex(Index bj) const
    {
      return IsColMajor ? blockOuterIndex(bj) : blockInnerIndex(bj);
    }

    inline Index blockOuterIndex(Index bj) const
    {
      return (m_blockSize == Dynamic) ? m_outerOffset[bj] : (bj * m_blockSize);
    }
    inline Index blockInnerIndex(Index bi) const
    {
      return (m_blockSize == Dynamic) ? m_innerOffset[bi] : (bi * m_blockSize);
    }

    // Not needed ???
    inline Index blockInnerSize(Index bi) const
    {
      return (m_blockSize == Dynamic) ? (m_innerOffset[bi+1] - m_innerOffset[bi]) : m_blockSize;
    }
    inline Index blockOuterSize(Index bj) const
    {
      return (m_blockSize == Dynamic) ? (m_outerOffset[bj+1]- m_outerOffset[bj]) : m_blockSize;
    }

    /**
      * \brief Browse the matrix by outer index
      */
    class InnerIterator; // Browse column by column

    /**
      * \brief Browse the matrix by block outer index
      */
    class BlockInnerIterator; // Browse block by block

    friend std::ostream & operator << (std::ostream & s, const BlockSparseMatrix& m)
    {
      for (StorageIndex j = 0; j < m.outerBlocks(); ++j)
      {
        BlockInnerIterator itb(m, j);
        for(; itb; ++itb)
        {
          s << "("<<itb.row() << ", " << itb.col() << ")\n";
          s << itb.value() <<"\n";
        }
      }
      s << std::endl;
      return s;
    }

    /**
      * \returns the starting position of the block \p id in the array of values
      */
    Index blockPtr(Index id) const
    {
      if(m_blockSize == Dynamic) return m_blockPtr[id];
      else return id * m_blockSize * m_blockSize;
      //return blockDynIdx(id, typename internal::conditional<(BlockSize==Dynamic), internal::true_type, internal::false_type>::type());
    }

  protected:
//    inline Index blockDynIdx(Index id, internal::true_type) const
//    {
//      return m_blockPtr[id];
//    }
//    inline Index blockDynIdx(Index id, internal::false_type) const
//    {
//      return id * BlockSize * BlockSize;
//    }

    // To be implemented
    // Insert a block at a particular location... need to make a room for that
    Map<BlockScalar> insert(Index brow, Index bcol);

    Index m_innerBSize; // Number of block rows
    Index m_outerBSize; // Number of block columns
    StorageIndex *m_innerOffset; // Starting index of each inner block (size m_innerBSize+1)
    StorageIndex *m_outerOffset; // Starting index of each outer block (size m_outerBSize+1)
    Index m_nonzerosblocks; // Total nonzeros blocks (lower than  m_innerBSize x m_outerBSize)
    Index m_nonzeros; // Total nonzeros elements
    Scalar *m_values; //Values stored block column after block column (size m_nonzeros)
    StorageIndex *m_blockPtr; // Pointer to the beginning of each block in m_values, size m_nonzeroblocks ... null for fixed-size blocks
    StorageIndex *m_indices; //Inner block indices, size m_nonzerosblocks ... OK
    StorageIndex *m_outerIndex; // Starting pointer of each block column in m_indices (size m_outerBSize)... OK
    Index m_blockSize; // Size of a block for fixed-size blocks, otherwise -1
};

template<typename _Scalar, int _BlockAtCompileTime, int _Options, typename _StorageIndex>
class BlockSparseMatrix<_Scalar, _BlockAtCompileTime, _Options, _StorageIndex>::BlockInnerIterator
{
  public:

    enum{
      Flags = _Options
    };

    BlockInnerIterator(const BlockSparseMatrix& mat, const Index outer)
    : m_mat(mat),m_outer(outer),
      m_id(mat.m_outerIndex[outer]),
      m_end(mat.m_outerIndex[outer+1])
    {
    }

    inline BlockInnerIterator& operator++() {m_id++; return *this; }

    inline const Map<const BlockScalar> value() const
    {
      return Map<const BlockScalar>(&(m_mat.m_values[m_mat.blockPtr(m_id)]),
          rows(),cols());
    }
    inline Map<BlockScalar> valueRef()
    {
      return Map<BlockScalar>(&(m_mat.m_values[m_mat.blockPtr(m_id)]),
          rows(),cols());
    }
    // Block inner index
    inline Index index() const {return m_mat.m_indices[m_id]; }
    inline Index outer() const { return m_outer; }
    // block row index
    inline Index row() const  {return index(); }
    // block column index
    inline Index col() const {return outer(); }
    // FIXME Number of rows in the current block
    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; }
    // Number of columns in the current block ...
    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;}
    inline operator bool() const { return (m_id < m_end); }

  protected:
    const BlockSparseMatrix<_Scalar, _BlockAtCompileTime, _Options, StorageIndex>& m_mat;
    const Index m_outer;
    Index m_id;
    Index m_end;
};

template<typename _Scalar, int _BlockAtCompileTime, int _Options, typename _StorageIndex>
class BlockSparseMatrix<_Scalar, _BlockAtCompileTime, _Options, _StorageIndex>::InnerIterator
{
  public:
    InnerIterator(const BlockSparseMatrix& mat, Index outer)
    : m_mat(mat),m_outerB(mat.outerToBlock(outer)),m_outer(outer),
      itb(mat, mat.outerToBlock(outer)),
      m_offset(outer - mat.blockOuterIndex(m_outerB))
     {
        if (itb)
        {
          m_id = m_mat.blockInnerIndex(itb.index());
          m_start = m_id;
          m_end = m_mat.blockInnerIndex(itb.index()+1);
        }
     }
    inline InnerIterator& operator++()
    {
      m_id++;
      if (m_id >= m_end)
      {
        ++itb;
        if (itb)
        {
          m_id = m_mat.blockInnerIndex(itb.index());
          m_start = m_id;
          m_end = m_mat.blockInnerIndex(itb.index()+1);
        }
      }
      return *this;
    }
    inline const Scalar& value() const
    {
      return itb.value().coeff(m_id - m_start, m_offset);
    }
    inline Scalar& valueRef()
    {
      return itb.valueRef().coeff(m_id - m_start, m_offset);
    }
    inline Index index() const { return m_id; }
    inline Index outer() const {return m_outer; }
    inline Index col() const {return outer(); }
    inline Index row() const { return index();}
    inline operator bool() const
    {
      return itb;
    }
  protected:
    const BlockSparseMatrix& m_mat;
    const Index m_outer;
    const Index m_outerB;
    BlockInnerIterator itb; // Iterator through the blocks
    const Index m_offset; // Position of this column in the block
    Index m_start; // starting inner index of this block
    Index m_id; // current inner index in the block
    Index m_end; // starting inner index of the next block

};
} // end namespace Eigen

#endif // EIGEN_SPARSEBLOCKMATRIX_H

quality91%

¤ Dauer der Verarbeitung: 0.4 Sekunden (vorverarbeitet) ¤

Wurzel

Suchen

Beweissystem der NASA

Beweissystem Isabelle

NIST Cobol Testsuite

Cephes Mathematical Library

Wiener Entwicklungsmethode

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung ist noch experimentell.