cart-elc

Source code for CART-ELC
git clone git://git.laack.co/cart-elc.git
Log | Files | Refs | README | LICENSE

SparseMatrixBase.h (17451B)


      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2008-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
      5 //
      6 // This Source Code Form is subject to the terms of the Mozilla
      7 // Public License v. 2.0. If a copy of the MPL was not distributed
      8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
      9 
     10 #ifndef EIGEN_SPARSEMATRIXBASE_H
     11 #define EIGEN_SPARSEMATRIXBASE_H
     12 
     13 namespace Eigen { 
     14 
     15 /** \ingroup SparseCore_Module
     16   *
     17   * \class SparseMatrixBase
     18   *
     19   * \brief Base class of any sparse matrices or sparse expressions
     20   *
     21   * \tparam Derived is the derived type, e.g. a sparse matrix type, or an expression, etc.
     22   *
     23   * This class can be extended with the help of the plugin mechanism described on the page
     24   * \ref TopicCustomizing_Plugins by defining the preprocessor symbol \c EIGEN_SPARSEMATRIXBASE_PLUGIN.
     25   */
     26 template<typename Derived> class SparseMatrixBase
     27   : public EigenBase<Derived>
     28 {
     29   public:
     30 
     31     typedef typename internal::traits<Derived>::Scalar Scalar;
     32     
     33     /** The numeric type of the expression' coefficients, e.g. float, double, int or std::complex<float>, etc.
     34       *
     35       * It is an alias for the Scalar type */
     36     typedef Scalar value_type;
     37     
     38     typedef typename internal::packet_traits<Scalar>::type PacketScalar;
     39     typedef typename internal::traits<Derived>::StorageKind StorageKind;
     40 
     41     /** The integer type used to \b store indices within a SparseMatrix.
     42       * For a \c SparseMatrix<Scalar,Options,IndexType> it an alias of the third template parameter \c IndexType. */
     43     typedef typename internal::traits<Derived>::StorageIndex StorageIndex;
     44 
     45     typedef typename internal::add_const_on_value_type_if_arithmetic<
     46                          typename internal::packet_traits<Scalar>::type
     47                      >::type PacketReturnType;
     48 
     49     typedef SparseMatrixBase StorageBaseType;
     50 
     51     typedef Matrix<StorageIndex,Dynamic,1> IndexVector;
     52     typedef Matrix<Scalar,Dynamic,1> ScalarVector;
     53     
     54     template<typename OtherDerived>
     55     Derived& operator=(const EigenBase<OtherDerived> &other);
     56 
     57     enum {
     58 
     59       RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime,
     60         /**< The number of rows at compile-time. This is just a copy of the value provided
     61           * by the \a Derived type. If a value is not known at compile-time,
     62           * it is set to the \a Dynamic constant.
     63           * \sa MatrixBase::rows(), MatrixBase::cols(), ColsAtCompileTime, SizeAtCompileTime */
     64 
     65       ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime,
     66         /**< The number of columns at compile-time. This is just a copy of the value provided
     67           * by the \a Derived type. If a value is not known at compile-time,
     68           * it is set to the \a Dynamic constant.
     69           * \sa MatrixBase::rows(), MatrixBase::cols(), RowsAtCompileTime, SizeAtCompileTime */
     70 
     71 
     72       SizeAtCompileTime = (internal::size_at_compile_time<internal::traits<Derived>::RowsAtCompileTime,
     73                                                    internal::traits<Derived>::ColsAtCompileTime>::ret),
     74         /**< This is equal to the number of coefficients, i.e. the number of
     75           * rows times the number of columns, or to \a Dynamic if this is not
     76           * known at compile-time. \sa RowsAtCompileTime, ColsAtCompileTime */
     77 
     78       MaxRowsAtCompileTime = RowsAtCompileTime,
     79       MaxColsAtCompileTime = ColsAtCompileTime,
     80 
     81       MaxSizeAtCompileTime = (internal::size_at_compile_time<MaxRowsAtCompileTime,
     82                                                       MaxColsAtCompileTime>::ret),
     83 
     84       IsVectorAtCompileTime = RowsAtCompileTime == 1 || ColsAtCompileTime == 1,
     85         /**< This is set to true if either the number of rows or the number of
     86           * columns is known at compile-time to be equal to 1. Indeed, in that case,
     87           * we are dealing with a column-vector (if there is only one column) or with
     88           * a row-vector (if there is only one row). */
     89 
     90       NumDimensions = int(MaxSizeAtCompileTime) == 1 ? 0 : bool(IsVectorAtCompileTime) ? 1 : 2,
     91         /**< This value is equal to Tensor::NumDimensions, i.e. 0 for scalars, 1 for vectors,
     92          * and 2 for matrices.
     93          */
     94 
     95       Flags = internal::traits<Derived>::Flags,
     96         /**< This stores expression \ref flags flags which may or may not be inherited by new expressions
     97           * constructed from this one. See the \ref flags "list of flags".
     98           */
     99 
    100       IsRowMajor = Flags&RowMajorBit ? 1 : 0,
    101       
    102       InnerSizeAtCompileTime = int(IsVectorAtCompileTime) ? int(SizeAtCompileTime)
    103                              : int(IsRowMajor) ? int(ColsAtCompileTime) : int(RowsAtCompileTime),
    104 
    105       #ifndef EIGEN_PARSED_BY_DOXYGEN
    106       _HasDirectAccess = (int(Flags)&DirectAccessBit) ? 1 : 0 // workaround sunCC
    107       #endif
    108     };
    109 
    110     /** \internal the return type of MatrixBase::adjoint() */
    111     typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
    112                         CwiseUnaryOp<internal::scalar_conjugate_op<Scalar>, Eigen::Transpose<const Derived> >,
    113                         Transpose<const Derived>
    114                      >::type AdjointReturnType;
    115     typedef Transpose<Derived> TransposeReturnType;
    116     typedef typename internal::add_const<Transpose<const Derived> >::type ConstTransposeReturnType;
    117 
    118     // FIXME storage order do not match evaluator storage order
    119     typedef SparseMatrix<Scalar, Flags&RowMajorBit ? RowMajor : ColMajor, StorageIndex> PlainObject;
    120 
    121 #ifndef EIGEN_PARSED_BY_DOXYGEN
    122     /** This is the "real scalar" type; if the \a Scalar type is already real numbers
    123       * (e.g. int, float or double) then \a RealScalar is just the same as \a Scalar. If
    124       * \a Scalar is \a std::complex<T> then RealScalar is \a T.
    125       *
    126       * \sa class NumTraits
    127       */
    128     typedef typename NumTraits<Scalar>::Real RealScalar;
    129 
    130     /** \internal the return type of coeff()
    131       */
    132     typedef typename internal::conditional<_HasDirectAccess, const Scalar&, Scalar>::type CoeffReturnType;
    133 
    134     /** \internal Represents a matrix with all coefficients equal to one another*/
    135     typedef CwiseNullaryOp<internal::scalar_constant_op<Scalar>,Matrix<Scalar,Dynamic,Dynamic> > ConstantReturnType;
    136 
    137     /** type of the equivalent dense matrix */
    138     typedef Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime> DenseMatrixType;
    139     /** type of the equivalent square matrix */
    140     typedef Matrix<Scalar,EIGEN_SIZE_MAX(RowsAtCompileTime,ColsAtCompileTime),
    141                           EIGEN_SIZE_MAX(RowsAtCompileTime,ColsAtCompileTime)> SquareMatrixType;
    142 
    143     inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
    144     inline Derived& derived() { return *static_cast<Derived*>(this); }
    145     inline Derived& const_cast_derived() const
    146     { return *static_cast<Derived*>(const_cast<SparseMatrixBase*>(this)); }
    147 
    148     typedef EigenBase<Derived> Base;
    149 
    150 #endif // not EIGEN_PARSED_BY_DOXYGEN
    151 
    152 #define EIGEN_CURRENT_STORAGE_BASE_CLASS Eigen::SparseMatrixBase
    153 #ifdef EIGEN_PARSED_BY_DOXYGEN
    154 #define EIGEN_DOC_UNARY_ADDONS(METHOD,OP)           /** <p>This method does not change the sparsity of \c *this: the OP is applied to explicitly stored coefficients only. \sa SparseCompressedBase::coeffs() </p> */
    155 #define EIGEN_DOC_BLOCK_ADDONS_NOT_INNER_PANEL      /** <p> \warning This method returns a read-only expression for any sparse matrices. \sa \ref TutorialSparse_SubMatrices "Sparse block operations" </p> */
    156 #define EIGEN_DOC_BLOCK_ADDONS_INNER_PANEL_IF(COND) /** <p> \warning This method returns a read-write expression for COND sparse matrices only. Otherwise, the returned expression is read-only. \sa \ref TutorialSparse_SubMatrices "Sparse block operations" </p> */
    157 #else
    158 #define EIGEN_DOC_UNARY_ADDONS(X,Y)
    159 #define EIGEN_DOC_BLOCK_ADDONS_NOT_INNER_PANEL
    160 #define EIGEN_DOC_BLOCK_ADDONS_INNER_PANEL_IF(COND)
    161 #endif
    162 #   include "../plugins/CommonCwiseUnaryOps.h"
    163 #   include "../plugins/CommonCwiseBinaryOps.h"
    164 #   include "../plugins/MatrixCwiseUnaryOps.h"
    165 #   include "../plugins/MatrixCwiseBinaryOps.h"
    166 #   include "../plugins/BlockMethods.h"
    167 #   ifdef EIGEN_SPARSEMATRIXBASE_PLUGIN
    168 #     include EIGEN_SPARSEMATRIXBASE_PLUGIN
    169 #   endif
    170 #undef EIGEN_CURRENT_STORAGE_BASE_CLASS
    171 #undef EIGEN_DOC_UNARY_ADDONS
    172 #undef EIGEN_DOC_BLOCK_ADDONS_NOT_INNER_PANEL
    173 #undef EIGEN_DOC_BLOCK_ADDONS_INNER_PANEL_IF
    174 
    175     /** \returns the number of rows. \sa cols() */
    176     inline Index rows() const { return derived().rows(); }
    177     /** \returns the number of columns. \sa rows() */
    178     inline Index cols() const { return derived().cols(); }
    179     /** \returns the number of coefficients, which is \a rows()*cols().
    180       * \sa rows(), cols(). */
    181     inline Index size() const { return rows() * cols(); }
    182     /** \returns true if either the number of rows or the number of columns is equal to 1.
    183       * In other words, this function returns
    184       * \code rows()==1 || cols()==1 \endcode
    185       * \sa rows(), cols(), IsVectorAtCompileTime. */
    186     inline bool isVector() const { return rows()==1 || cols()==1; }
    187     /** \returns the size of the storage major dimension,
    188       * i.e., the number of columns for a columns major matrix, and the number of rows otherwise */
    189     Index outerSize() const { return (int(Flags)&RowMajorBit) ? this->rows() : this->cols(); }
    190     /** \returns the size of the inner dimension according to the storage order,
    191       * i.e., the number of rows for a columns major matrix, and the number of cols otherwise */
    192     Index innerSize() const { return (int(Flags)&RowMajorBit) ? this->cols() : this->rows(); }
    193 
    194     bool isRValue() const { return m_isRValue; }
    195     Derived& markAsRValue() { m_isRValue = true; return derived(); }
    196 
    197     SparseMatrixBase() : m_isRValue(false) { /* TODO check flags */ }
    198 
    199     
    200     template<typename OtherDerived>
    201     Derived& operator=(const ReturnByValue<OtherDerived>& other);
    202 
    203     template<typename OtherDerived>
    204     inline Derived& operator=(const SparseMatrixBase<OtherDerived>& other);
    205 
    206     inline Derived& operator=(const Derived& other);
    207 
    208   protected:
    209 
    210     template<typename OtherDerived>
    211     inline Derived& assign(const OtherDerived& other);
    212 
    213     template<typename OtherDerived>
    214     inline void assignGeneric(const OtherDerived& other);
    215 
    216   public:
    217 
    218     friend std::ostream & operator << (std::ostream & s, const SparseMatrixBase& m)
    219     {
    220       typedef typename Derived::Nested Nested;
    221       typedef typename internal::remove_all<Nested>::type NestedCleaned;
    222 
    223       if (Flags&RowMajorBit)
    224       {
    225         Nested nm(m.derived());
    226         internal::evaluator<NestedCleaned> thisEval(nm);
    227         for (Index row=0; row<nm.outerSize(); ++row)
    228         {
    229           Index col = 0;
    230           for (typename internal::evaluator<NestedCleaned>::InnerIterator it(thisEval, row); it; ++it)
    231           {
    232             for ( ; col<it.index(); ++col)
    233               s << "0 ";
    234             s << it.value() << " ";
    235             ++col;
    236           }
    237           for ( ; col<m.cols(); ++col)
    238             s << "0 ";
    239           s << std::endl;
    240         }
    241       }
    242       else
    243       {
    244         Nested nm(m.derived());
    245         internal::evaluator<NestedCleaned> thisEval(nm);
    246         if (m.cols() == 1) {
    247           Index row = 0;
    248           for (typename internal::evaluator<NestedCleaned>::InnerIterator it(thisEval, 0); it; ++it)
    249           {
    250             for ( ; row<it.index(); ++row)
    251               s << "0" << std::endl;
    252             s << it.value() << std::endl;
    253             ++row;
    254           }
    255           for ( ; row<m.rows(); ++row)
    256             s << "0" << std::endl;
    257         }
    258         else
    259         {
    260           SparseMatrix<Scalar, RowMajorBit, StorageIndex> trans = m;
    261           s << static_cast<const SparseMatrixBase<SparseMatrix<Scalar, RowMajorBit, StorageIndex> >&>(trans);
    262         }
    263       }
    264       return s;
    265     }
    266 
    267     template<typename OtherDerived>
    268     Derived& operator+=(const SparseMatrixBase<OtherDerived>& other);
    269     template<typename OtherDerived>
    270     Derived& operator-=(const SparseMatrixBase<OtherDerived>& other);
    271     
    272     template<typename OtherDerived>
    273     Derived& operator+=(const DiagonalBase<OtherDerived>& other);
    274     template<typename OtherDerived>
    275     Derived& operator-=(const DiagonalBase<OtherDerived>& other);
    276 
    277     template<typename OtherDerived>
    278     Derived& operator+=(const EigenBase<OtherDerived> &other);
    279     template<typename OtherDerived>
    280     Derived& operator-=(const EigenBase<OtherDerived> &other);
    281 
    282     Derived& operator*=(const Scalar& other);
    283     Derived& operator/=(const Scalar& other);
    284 
    285     template<typename OtherDerived> struct CwiseProductDenseReturnType {
    286       typedef CwiseBinaryOp<internal::scalar_product_op<typename ScalarBinaryOpTraits<
    287                                                           typename internal::traits<Derived>::Scalar,
    288                                                           typename internal::traits<OtherDerived>::Scalar
    289                                                         >::ReturnType>,
    290                             const Derived,
    291                             const OtherDerived
    292                           > Type;
    293     };
    294 
    295     template<typename OtherDerived>
    296     EIGEN_STRONG_INLINE const typename CwiseProductDenseReturnType<OtherDerived>::Type
    297     cwiseProduct(const MatrixBase<OtherDerived> &other) const;
    298 
    299     // sparse * diagonal
    300     template<typename OtherDerived>
    301     const Product<Derived,OtherDerived>
    302     operator*(const DiagonalBase<OtherDerived> &other) const
    303     { return Product<Derived,OtherDerived>(derived(), other.derived()); }
    304 
    305     // diagonal * sparse
    306     template<typename OtherDerived> friend
    307     const Product<OtherDerived,Derived>
    308     operator*(const DiagonalBase<OtherDerived> &lhs, const SparseMatrixBase& rhs)
    309     { return Product<OtherDerived,Derived>(lhs.derived(), rhs.derived()); }
    310     
    311     // sparse * sparse
    312     template<typename OtherDerived>
    313     const Product<Derived,OtherDerived,AliasFreeProduct>
    314     operator*(const SparseMatrixBase<OtherDerived> &other) const;
    315     
    316     // sparse * dense
    317     template<typename OtherDerived>
    318     const Product<Derived,OtherDerived>
    319     operator*(const MatrixBase<OtherDerived> &other) const
    320     { return Product<Derived,OtherDerived>(derived(), other.derived()); }
    321     
    322     // dense * sparse
    323     template<typename OtherDerived> friend
    324     const Product<OtherDerived,Derived>
    325     operator*(const MatrixBase<OtherDerived> &lhs, const SparseMatrixBase& rhs)
    326     { return Product<OtherDerived,Derived>(lhs.derived(), rhs.derived()); }
    327     
    328      /** \returns an expression of P H P^-1 where H is the matrix represented by \c *this */
    329     SparseSymmetricPermutationProduct<Derived,Upper|Lower> twistedBy(const PermutationMatrix<Dynamic,Dynamic,StorageIndex>& perm) const
    330     {
    331       return SparseSymmetricPermutationProduct<Derived,Upper|Lower>(derived(), perm);
    332     }
    333 
    334     template<typename OtherDerived>
    335     Derived& operator*=(const SparseMatrixBase<OtherDerived>& other);
    336 
    337     template<int Mode>
    338     inline const TriangularView<const Derived, Mode> triangularView() const;
    339     
    340     template<unsigned int UpLo> struct SelfAdjointViewReturnType { typedef SparseSelfAdjointView<Derived, UpLo> Type; };
    341     template<unsigned int UpLo> struct ConstSelfAdjointViewReturnType { typedef const SparseSelfAdjointView<const Derived, UpLo> Type; };
    342 
    343     template<unsigned int UpLo> inline 
    344     typename ConstSelfAdjointViewReturnType<UpLo>::Type selfadjointView() const;
    345     template<unsigned int UpLo> inline
    346     typename SelfAdjointViewReturnType<UpLo>::Type selfadjointView();
    347 
    348     template<typename OtherDerived> Scalar dot(const MatrixBase<OtherDerived>& other) const;
    349     template<typename OtherDerived> Scalar dot(const SparseMatrixBase<OtherDerived>& other) const;
    350     RealScalar squaredNorm() const;
    351     RealScalar norm()  const;
    352     RealScalar blueNorm() const;
    353 
    354     TransposeReturnType transpose() { return TransposeReturnType(derived()); }
    355     const ConstTransposeReturnType transpose() const { return ConstTransposeReturnType(derived()); }
    356     const AdjointReturnType adjoint() const { return AdjointReturnType(transpose()); }
    357 
    358     DenseMatrixType toDense() const
    359     {
    360       return DenseMatrixType(derived());
    361     }
    362 
    363     template<typename OtherDerived>
    364     bool isApprox(const SparseMatrixBase<OtherDerived>& other,
    365                   const RealScalar& prec = NumTraits<Scalar>::dummy_precision()) const;
    366 
    367     template<typename OtherDerived>
    368     bool isApprox(const MatrixBase<OtherDerived>& other,
    369                   const RealScalar& prec = NumTraits<Scalar>::dummy_precision()) const
    370     { return toDense().isApprox(other,prec); }
    371 
    372     /** \returns the matrix or vector obtained by evaluating this expression.
    373       *
    374       * Notice that in the case of a plain matrix or vector (not an expression) this function just returns
    375       * a const reference, in order to avoid a useless copy.
    376       */
    377     inline const typename internal::eval<Derived>::type eval() const
    378     { return typename internal::eval<Derived>::type(derived()); }
    379 
    380     Scalar sum() const;
    381     
    382     inline const SparseView<Derived>
    383     pruned(const Scalar& reference = Scalar(0), const RealScalar& epsilon = NumTraits<Scalar>::dummy_precision()) const;
    384 
    385   protected:
    386 
    387     bool m_isRValue;
    388 
    389     static inline StorageIndex convert_index(const Index idx) {
    390       return internal::convert_index<StorageIndex>(idx);
    391     }
    392   private:
    393     template<typename Dest> void evalTo(Dest &) const;
    394 };
    395 
    396 } // end namespace Eigen
    397 
    398 #endif // EIGEN_SPARSEMATRIXBASE_H