// This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2008-2009 Gael Guennebaud // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public // License as published by the Free Software Foundation; either // version 3 of the License, or (at your option) any later version. // // Alternatively, you can redistribute it and/or // modify it under the terms of the GNU General Public License as // published by the Free Software Foundation; either version 2 of // the License, or (at your option) any later version. // // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the // GNU General Public License for more details. // // You should have received a copy of the GNU Lesser General Public // License and a copy of the GNU General Public License along with // Eigen. If not, see . #ifndef EIGEN_SPARSEMATRIXBASE_H #define EIGEN_SPARSEMATRIXBASE_H /** \ingroup Sparse_Module * * \class SparseMatrixBase * * \brief Base class of any sparse matrices or sparse expressions * * \param Derived * * * */ template class SparseMatrixBase : public EigenBase { public: typedef typename ei_traits::Scalar Scalar; typedef typename ei_packet_traits::type PacketScalar; typedef typename ei_traits::StorageKind StorageKind; typedef typename ei_traits::Index Index; typedef SparseMatrixBase StorageBaseType; enum { RowsAtCompileTime = ei_traits::RowsAtCompileTime, /**< The number of rows at compile-time. This is just a copy of the value provided * by the \a Derived type. If a value is not known at compile-time, * it is set to the \a Dynamic constant. * \sa MatrixBase::rows(), MatrixBase::cols(), ColsAtCompileTime, SizeAtCompileTime */ ColsAtCompileTime = ei_traits::ColsAtCompileTime, /**< The number of columns at compile-time. This is just a copy of the value provided * by the \a Derived type. If a value is not known at compile-time, * it is set to the \a Dynamic constant. * \sa MatrixBase::rows(), MatrixBase::cols(), RowsAtCompileTime, SizeAtCompileTime */ SizeAtCompileTime = (ei_size_at_compile_time::RowsAtCompileTime, ei_traits::ColsAtCompileTime>::ret), /**< This is equal to the number of coefficients, i.e. the number of * rows times the number of columns, or to \a Dynamic if this is not * known at compile-time. \sa RowsAtCompileTime, ColsAtCompileTime */ MaxRowsAtCompileTime = RowsAtCompileTime, MaxColsAtCompileTime = ColsAtCompileTime, MaxSizeAtCompileTime = (ei_size_at_compile_time::ret), IsVectorAtCompileTime = RowsAtCompileTime == 1 || ColsAtCompileTime == 1, /**< This is set to true if either the number of rows or the number of * columns is known at compile-time to be equal to 1. Indeed, in that case, * we are dealing with a column-vector (if there is only one column) or with * a row-vector (if there is only one row). */ Flags = ei_traits::Flags, /**< This stores expression \ref flags flags which may or may not be inherited by new expressions * constructed from this one. See the \ref flags "list of flags". */ CoeffReadCost = ei_traits::CoeffReadCost, /**< This is a rough measure of how expensive it is to read one coefficient from * this expression. */ IsRowMajor = Flags&RowMajorBit ? 1 : 0, #ifndef EIGEN_PARSED_BY_DOXYGEN _HasDirectAccess = (int(Flags)&DirectAccessBit) ? 1 : 0 // workaround sunCC #endif }; /* \internal the return type of MatrixBase::conjugate() */ // typedef typename ei_meta_if::IsComplex, // const SparseCwiseUnaryOp, Derived>, // const Derived& // >::ret ConjugateReturnType; /* \internal the return type of MatrixBase::real() */ // typedef SparseCwiseUnaryOp, Derived> RealReturnType; /* \internal the return type of MatrixBase::imag() */ // typedef SparseCwiseUnaryOp, Derived> ImagReturnType; /** \internal the return type of MatrixBase::adjoint() */ typedef typename ei_meta_if::IsComplex, CwiseUnaryOp, Eigen::Transpose >, Transpose >::ret AdjointReturnType; typedef SparseMatrix PlainObject; #define EIGEN_CURRENT_STORAGE_BASE_CLASS Eigen::SparseMatrixBase #include "../plugins/CommonCwiseUnaryOps.h" #include "../plugins/CommonCwiseBinaryOps.h" #include "../plugins/MatrixCwiseUnaryOps.h" #include "../plugins/MatrixCwiseBinaryOps.h" #undef EIGEN_CURRENT_STORAGE_BASE_CLASS #ifndef EIGEN_PARSED_BY_DOXYGEN /** This is the "real scalar" type; if the \a Scalar type is already real numbers * (e.g. int, float or double) then \a RealScalar is just the same as \a Scalar. If * \a Scalar is \a std::complex then RealScalar is \a T. * * \sa class NumTraits */ typedef typename NumTraits::Real RealScalar; /** \internal the return type of coeff() */ typedef typename ei_meta_if<_HasDirectAccess, const Scalar&, Scalar>::ret CoeffReturnType; /** \internal Represents a matrix with all coefficients equal to one another*/ typedef CwiseNullaryOp,Matrix > ConstantReturnType; /** type of the equivalent square matrix */ typedef Matrix SquareMatrixType; inline const Derived& derived() const { return *static_cast(this); } inline Derived& derived() { return *static_cast(this); } inline Derived& const_cast_derived() const { return *static_cast(const_cast(this)); } #endif // not EIGEN_PARSED_BY_DOXYGEN /** \returns the number of rows. \sa cols(), RowsAtCompileTime */ inline Index rows() const { return derived().rows(); } /** \returns the number of columns. \sa rows(), ColsAtCompileTime*/ inline Index cols() const { return derived().cols(); } /** \returns the number of coefficients, which is \a rows()*cols(). * \sa rows(), cols(), SizeAtCompileTime. */ inline Index size() const { return rows() * cols(); } /** \returns the number of nonzero coefficients which is in practice the number * of stored coefficients. */ inline Index nonZeros() const { return derived().nonZeros(); } /** \returns true if either the number of rows or the number of columns is equal to 1. * In other words, this function returns * \code rows()==1 || cols()==1 \endcode * \sa rows(), cols(), IsVectorAtCompileTime. */ inline bool isVector() const { return rows()==1 || cols()==1; } /** \returns the size of the storage major dimension, * i.e., the number of columns for a columns major matrix, and the number of rows otherwise */ Index outerSize() const { return (int(Flags)&RowMajorBit) ? this->rows() : this->cols(); } /** \returns the size of the inner dimension according to the storage order, * i.e., the number of rows for a columns major matrix, and the number of cols otherwise */ Index innerSize() const { return (int(Flags)&RowMajorBit) ? this->cols() : this->rows(); } bool isRValue() const { return m_isRValue; } Derived& markAsRValue() { m_isRValue = true; return derived(); } SparseMatrixBase() : m_isRValue(false) { /* TODO check flags */ } inline Derived& operator=(const Derived& other) { // std::cout << "Derived& operator=(const Derived& other)\n"; // if (other.isRValue()) // derived().swap(other.const_cast_derived()); // else this->operator=(other); return derived(); } template inline void assignGeneric(const OtherDerived& other) { // std::cout << "Derived& operator=(const MatrixBase& other)\n"; //const bool transpose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit); ei_assert(( ((ei_traits::SupportedAccessPatterns&OuterRandomAccessPattern)==OuterRandomAccessPattern) || (!((Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit)))) && "the transpose operation is supposed to be handled in SparseMatrix::operator="); enum { Flip = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit) }; const Index outerSize = other.outerSize(); //typedef typename ei_meta_if, Derived>::ret TempType; // thanks to shallow copies, we always eval to a tempary Derived temp(other.rows(), other.cols()); temp.reserve(std::max(this->rows(),this->cols())*2); for (Index j=0; j inline Derived& operator=(const SparseMatrixBase& other) { // std::cout << typeid(OtherDerived).name() << "\n"; // std::cout << Flags << " " << OtherDerived::Flags << "\n"; const bool transpose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit); // std::cout << "eval transpose = " << transpose << "\n"; const Index outerSize = (int(OtherDerived::Flags) & RowMajorBit) ? other.rows() : other.cols(); if ((!transpose) && other.isRValue()) { // eval without temporary derived().resize(other.rows(), other.cols()); derived().setZero(); derived().reserve(std::max(this->rows(),this->cols())*2); for (Index j=0; j inline Derived& operator=(const SparseSparseProduct& product); template inline void _experimentalNewProduct(const Lhs& lhs, const Rhs& rhs); friend std::ostream & operator << (std::ostream & s, const SparseMatrixBase& m) { if (Flags&RowMajorBit) { for (Index row=0; row trans = m.derived(); s << trans; } } return s; } // const SparseCwiseUnaryOp::Scalar>,Derived> operator-() const; // template // const CwiseBinaryOp::Scalar>, Derived, OtherDerived> // operator+(const SparseMatrixBase &other) const; // template // const CwiseBinaryOp::Scalar>, Derived, OtherDerived> // operator-(const SparseMatrixBase &other) const; template Derived& operator+=(const SparseMatrixBase& other); template Derived& operator-=(const SparseMatrixBase& other); // template // Derived& operator+=(const Flagged, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit>& other); Derived& operator*=(const Scalar& other); Derived& operator/=(const Scalar& other); #define EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE \ CwiseBinaryOp< \ ei_scalar_product_op< \ typename ei_scalar_product_traits< \ typename ei_traits::Scalar, \ typename ei_traits::Scalar \ >::ReturnType \ >, \ Derived, \ OtherDerived \ > template EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE cwiseProduct(const MatrixBase &other) const; // const SparseCwiseUnaryOp::Scalar>, Derived> // operator*(const Scalar& scalar) const; // const SparseCwiseUnaryOp::Scalar>, Derived> // operator/(const Scalar& scalar) const; // inline friend const SparseCwiseUnaryOp::Scalar>, Derived> // operator*(const Scalar& scalar, const SparseMatrixBase& matrix) // { return matrix*scalar; } // sparse * sparse template const typename SparseSparseProductReturnType::Type operator*(const SparseMatrixBase &other) const; // sparse * diagonal template const SparseDiagonalProduct operator*(const DiagonalBase &other) const; // diagonal * sparse template friend const SparseDiagonalProduct operator*(const DiagonalBase &lhs, const SparseMatrixBase& rhs) { return SparseDiagonalProduct(lhs.derived(), rhs.derived()); } /** dense * sparse (return a dense object unless it is an outer product) */ template friend const typename DenseSparseProductReturnType::Type operator*(const MatrixBase& lhs, const Derived& rhs) { return typename DenseSparseProductReturnType::Type(lhs.derived(),rhs); } /** sparse * dense (returns a dense object unless it is an outer product) */ template const typename SparseDenseProductReturnType::Type operator*(const MatrixBase &other) const; template Derived& operator*=(const SparseMatrixBase& other); #ifdef EIGEN2_SUPPORT // deprecated template typename ei_plain_matrix_type_column_major::type solveTriangular(const MatrixBase& other) const; // deprecated template void solveTriangularInPlace(MatrixBase& other) const; // template // void solveTriangularInPlace(SparseMatrixBase& other) const; #endif // EIGEN2_SUPPORT template inline const SparseTriangularView triangularView() const; template inline const SparseSelfAdjointView selfadjointView() const; template inline SparseSelfAdjointView selfadjointView(); template Scalar dot(const MatrixBase& other) const; template Scalar dot(const SparseMatrixBase& other) const; RealScalar squaredNorm() const; RealScalar norm() const; // const PlainObject normalized() const; // void normalize(); Transpose transpose() { return derived(); } const Transpose transpose() const { return derived(); } // void transposeInPlace(); const AdjointReturnType adjoint() const { return transpose(); } // sub-vector SparseInnerVectorSet row(Index i); const SparseInnerVectorSet row(Index i) const; SparseInnerVectorSet col(Index j); const SparseInnerVectorSet col(Index j) const; SparseInnerVectorSet innerVector(Index outer); const SparseInnerVectorSet innerVector(Index outer) const; // set of sub-vectors SparseInnerVectorSet subrows(Index start, Index size); const SparseInnerVectorSet subrows(Index start, Index size) const; SparseInnerVectorSet subcols(Index start, Index size); const SparseInnerVectorSet subcols(Index start, Index size) const; SparseInnerVectorSet innerVectors(Index outerStart, Index outerSize); const SparseInnerVectorSet innerVectors(Index outerStart, Index outerSize) const; // typename BlockReturnType::Type block(int startRow, int startCol, int blockRows, int blockCols); // const typename BlockReturnType::Type // block(int startRow, int startCol, int blockRows, int blockCols) const; // // typename BlockReturnType::SubVectorType segment(int start, int size); // const typename BlockReturnType::SubVectorType segment(int start, int size) const; // // typename BlockReturnType::SubVectorType start(int size); // const typename BlockReturnType::SubVectorType start(int size) const; // // typename BlockReturnType::SubVectorType end(int size); // const typename BlockReturnType::SubVectorType end(int size) const; // // template // typename BlockReturnType::Type block(int startRow, int startCol); // template // const typename BlockReturnType::Type block(int startRow, int startCol) const; // template typename BlockReturnType::SubVectorType start(void); // template const typename BlockReturnType::SubVectorType start() const; // template typename BlockReturnType::SubVectorType end(); // template const typename BlockReturnType::SubVectorType end() const; // template typename BlockReturnType::SubVectorType segment(int start); // template const typename BlockReturnType::SubVectorType segment(int start) const; // Diagonal diagonal(); // const Diagonal diagonal() const; // template Part part(); // template const Part part() const; // static const ConstantReturnType Constant(int rows, int cols, const Scalar& value); // static const ConstantReturnType Constant(int size, const Scalar& value); // static const ConstantReturnType Constant(const Scalar& value); // template // static const CwiseNullaryOp NullaryExpr(int rows, int cols, const CustomNullaryOp& func); // template // static const CwiseNullaryOp NullaryExpr(int size, const CustomNullaryOp& func); // template // static const CwiseNullaryOp NullaryExpr(const CustomNullaryOp& func); // static const ConstantReturnType Zero(int rows, int cols); // static const ConstantReturnType Zero(int size); // static const ConstantReturnType Zero(); // static const ConstantReturnType Ones(int rows, int cols); // static const ConstantReturnType Ones(int size); // static const ConstantReturnType Ones(); // static const IdentityReturnType Identity(); // static const IdentityReturnType Identity(int rows, int cols); // static const BasisReturnType Unit(int size, int i); // static const BasisReturnType Unit(int i); // static const BasisReturnType UnitX(); // static const BasisReturnType UnitY(); // static const BasisReturnType UnitZ(); // static const BasisReturnType UnitW(); // const DiagonalMatrix asDiagonal() const; // Derived& setConstant(const Scalar& value); // Derived& setZero(); // Derived& setOnes(); // Derived& setRandom(); // Derived& setIdentity(); /** \internal use operator= */ template void evalTo(MatrixBase& dst) const { dst.setZero(); for (Index j=0; j toDense() const { return derived(); } template bool isApprox(const SparseMatrixBase& other, RealScalar prec = NumTraits::dummy_precision()) const { return toDense().isApprox(other.toDense(),prec); } template bool isApprox(const MatrixBase& other, RealScalar prec = NumTraits::dummy_precision()) const { return toDense().isApprox(other,prec); } // bool isMuchSmallerThan(const RealScalar& other, // RealScalar prec = NumTraits::dummy_precision()) const; // template // bool isMuchSmallerThan(const MatrixBase& other, // RealScalar prec = NumTraits::dummy_precision()) const; // bool isApproxToConstant(const Scalar& value, RealScalar prec = NumTraits::dummy_precision()) const; // bool isZero(RealScalar prec = NumTraits::dummy_precision()) const; // bool isOnes(RealScalar prec = NumTraits::dummy_precision()) const; // bool isIdentity(RealScalar prec = NumTraits::dummy_precision()) const; // bool isDiagonal(RealScalar prec = NumTraits::dummy_precision()) const; // bool isUpper(RealScalar prec = NumTraits::dummy_precision()) const; // bool isLower(RealScalar prec = NumTraits::dummy_precision()) const; // template // bool isOrthogonal(const MatrixBase& other, // RealScalar prec = NumTraits::dummy_precision()) const; // bool isUnitary(RealScalar prec = NumTraits::dummy_precision()) const; // template // inline bool operator==(const MatrixBase& other) const // { return (cwise() == other).all(); } // template // inline bool operator!=(const MatrixBase& other) const // { return (cwise() != other).any(); } // template // const SparseCwiseUnaryOp::Scalar, NewType>, Derived> cast() const; /** \returns the matrix or vector obtained by evaluating this expression. * * Notice that in the case of a plain matrix or vector (not an expression) this function just returns * a const reference, in order to avoid a useless copy. */ inline const typename ei_eval::type eval() const { return typename ei_eval::type(derived()); } // template // void swap(MatrixBase EIGEN_REF_TO_TEMPORARY other); // template // const SparseFlagged marked() const; // const Flagged lazy() const; /** \returns number of elements to skip to pass from one row (resp. column) to another * for a row-major (resp. column-major) matrix. * Combined with coeffRef() and the \ref flags flags, it allows a direct access to the data * of the underlying matrix. */ // inline int stride(void) const { return derived().stride(); } // FIXME // ConjugateReturnType conjugate() const; // const RealReturnType real() const; // const ImagReturnType imag() const; // template // const SparseCwiseUnaryOp unaryExpr(const CustomUnaryOp& func = CustomUnaryOp()) const; // template // const CwiseBinaryOp // binaryExpr(const MatrixBase &other, const CustomBinaryOp& func = CustomBinaryOp()) const; Scalar sum() const; // Scalar trace() const; // typename ei_traits::Scalar minCoeff() const; // typename ei_traits::Scalar maxCoeff() const; // typename ei_traits::Scalar minCoeff(int* row, int* col = 0) const; // typename ei_traits::Scalar maxCoeff(int* row, int* col = 0) const; // template // typename ei_result_of::Scalar)>::type // redux(const BinaryOp& func) const; // template // void visit(Visitor& func) const; // const SparseCwise cwise() const; // SparseCwise cwise(); // inline const WithFormat format(const IOFormat& fmt) const; /////////// Array module /////////// /* bool all(void) const; bool any(void) const; const VectorwiseOp rowwise() const; const VectorwiseOp colwise() const; static const CwiseNullaryOp,Derived> Random(int rows, int cols); static const CwiseNullaryOp,Derived> Random(int size); static const CwiseNullaryOp,Derived> Random(); template const Select select(const MatrixBase& thenMatrix, const MatrixBase& elseMatrix) const; template inline const Select select(const MatrixBase& thenMatrix, typename ThenDerived::Scalar elseScalar) const; template inline const Select select(typename ElseDerived::Scalar thenScalar, const MatrixBase& elseMatrix) const; template RealScalar lpNorm() const; */ // template // Scalar dot(const MatrixBase& other) const // { // EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) // EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived) // EIGEN_STATIC_ASSERT((ei_is_same_type::ret), // YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) // // ei_assert(derived().size() == other.size()); // // short version, but the assembly looks more complicated because // // of the CwiseBinaryOp iterator complexity // // return res = (derived().cwise() * other.derived().conjugate()).sum(); // // // optimized, generic version // typename Derived::InnerIterator i(derived(),0); // typename OtherDerived::InnerIterator j(other.derived(),0); // Scalar res = 0; // while (i && j) // { // if (i.index()==j.index()) // { // // std::cerr << i.value() << " * " << j.value() << "\n"; // res += i.value() * ei_conj(j.value()); // ++i; ++j; // } // else if (i.index()