SparseMatrixBase.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_SPARSEMATRIXBASE_H
00026 #define EIGEN_SPARSEMATRIXBASE_H
00027 
00028 namespace Eigen { 
00029 
00041 template<typename Derived> class SparseMatrixBase : public EigenBase<Derived>
00042 {
00043   public:
00044 
00045     typedef typename internal::traits<Derived>::Scalar Scalar;
00046     typedef typename internal::packet_traits<Scalar>::type PacketScalar;
00047     typedef typename internal::traits<Derived>::StorageKind StorageKind;
00048     typedef typename internal::traits<Derived>::Index Index;
00049     typedef typename internal::add_const_on_value_type_if_arithmetic<
00050                          typename internal::packet_traits<Scalar>::type
00051                      >::type PacketReturnType;
00052 
00053     typedef SparseMatrixBase StorageBaseType;
00054     typedef EigenBase<Derived> Base;
00055     
00056     template<typename OtherDerived>
00057     Derived& operator=(const EigenBase<OtherDerived> &other)
00058     {
00059       other.derived().evalTo(derived());
00060       return derived();
00061     }
00062 
00063     enum {
00064 
00065       RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime,
00071       ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime,
00078       SizeAtCompileTime = (internal::size_at_compile_time<internal::traits<Derived>::RowsAtCompileTime,
00079                                                    internal::traits<Derived>::ColsAtCompileTime>::ret),
00084       MaxRowsAtCompileTime = RowsAtCompileTime,
00085       MaxColsAtCompileTime = ColsAtCompileTime,
00086 
00087       MaxSizeAtCompileTime = (internal::size_at_compile_time<MaxRowsAtCompileTime,
00088                                                       MaxColsAtCompileTime>::ret),
00089 
00090       IsVectorAtCompileTime = RowsAtCompileTime == 1 || ColsAtCompileTime == 1,
00096       Flags = internal::traits<Derived>::Flags,
00101       CoeffReadCost = internal::traits<Derived>::CoeffReadCost,
00106       IsRowMajor = Flags&RowMajorBit ? 1 : 0,
00107 
00108       #ifndef EIGEN_PARSED_BY_DOXYGEN
00109       _HasDirectAccess = (int(Flags)&DirectAccessBit) ? 1 : 0 // workaround sunCC
00110       #endif
00111     };
00112 
00114     typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
00115                         CwiseUnaryOp<internal::scalar_conjugate_op<Scalar>, Eigen::Transpose<const Derived> >,
00116                         Transpose<const Derived>
00117                      >::type AdjointReturnType;
00118 
00119 
00120     typedef SparseMatrix<Scalar, Flags&RowMajorBit ? RowMajor : ColMajor> PlainObject;
00121 
00122 
00123 #ifndef EIGEN_PARSED_BY_DOXYGEN
00124 
00130     typedef typename NumTraits<Scalar>::Real RealScalar;
00131 
00134     typedef typename internal::conditional<_HasDirectAccess, const Scalar&, Scalar>::type CoeffReturnType;
00135 
00137     typedef CwiseNullaryOp<internal::scalar_constant_op<Scalar>,Matrix<Scalar,Dynamic,Dynamic> > ConstantReturnType;
00138 
00140     typedef Matrix<Scalar,EIGEN_SIZE_MAX(RowsAtCompileTime,ColsAtCompileTime),
00141                           EIGEN_SIZE_MAX(RowsAtCompileTime,ColsAtCompileTime)> SquareMatrixType;
00142 
00143     inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
00144     inline Derived& derived() { return *static_cast<Derived*>(this); }
00145     inline Derived& const_cast_derived() const
00146     { return *static_cast<Derived*>(const_cast<SparseMatrixBase*>(this)); }
00147 #endif // not EIGEN_PARSED_BY_DOXYGEN
00148 
00149 #define EIGEN_CURRENT_STORAGE_BASE_CLASS Eigen::SparseMatrixBase
00150 #   include "../plugins/CommonCwiseUnaryOps.h"
00151 #   include "../plugins/CommonCwiseBinaryOps.h"
00152 #   include "../plugins/MatrixCwiseUnaryOps.h"
00153 #   include "../plugins/MatrixCwiseBinaryOps.h"
00154 #   ifdef EIGEN_SPARSEMATRIXBASE_PLUGIN
00155 #     include EIGEN_SPARSEMATRIXBASE_PLUGIN
00156 #   endif
00157 #   undef EIGEN_CURRENT_STORAGE_BASE_CLASS
00158 #undef EIGEN_CURRENT_STORAGE_BASE_CLASS
00159 
00160 
00162     inline Index rows() const { return derived().rows(); }
00164     inline Index cols() const { return derived().cols(); }
00167     inline Index size() const { return rows() * cols(); }
00170     inline Index nonZeros() const { return derived().nonZeros(); }
00175     inline bool isVector() const { return rows()==1 || cols()==1; }
00178     Index outerSize() const { return (int(Flags)&RowMajorBit) ? this->rows() : this->cols(); }
00181     Index innerSize() const { return (int(Flags)&RowMajorBit) ? this->cols() : this->rows(); }
00182 
00183     bool isRValue() const { return m_isRValue; }
00184     Derived& markAsRValue() { m_isRValue = true; return derived(); }
00185 
00186     SparseMatrixBase() : m_isRValue(false) { /* TODO check flags */ }
00187 
00188     
00189     template<typename OtherDerived>
00190     Derived& operator=(const ReturnByValue<OtherDerived>& other)
00191     {
00192       other.evalTo(derived());
00193       return derived();
00194     }
00195 
00196 
00197     template<typename OtherDerived>
00198     inline Derived& operator=(const SparseMatrixBase<OtherDerived>& other)
00199     {
00200       return assign(other.derived());
00201     }
00202 
00203     inline Derived& operator=(const Derived& other)
00204     {
00205 //       if (other.isRValue())
00206 //         derived().swap(other.const_cast_derived());
00207 //       else
00208       return assign(other.derived());
00209     }
00210 
00211   protected:
00212 
00213     template<typename OtherDerived>
00214     inline Derived& assign(const OtherDerived& other)
00215     {
00216       const bool transpose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit);
00217       const Index outerSize = (int(OtherDerived::Flags) & RowMajorBit) ? other.rows() : other.cols();
00218       if ((!transpose) && other.isRValue())
00219       {
00220         // eval without temporary
00221         derived().resize(other.rows(), other.cols());
00222         derived().setZero();
00223         derived().reserve((std::max)(this->rows(),this->cols())*2);
00224         for (Index j=0; j<outerSize; ++j)
00225         {
00226           derived().startVec(j);
00227           for (typename OtherDerived::InnerIterator it(other, j); it; ++it)
00228           {
00229             Scalar v = it.value();
00230             derived().insertBackByOuterInner(j,it.index()) = v;
00231           }
00232         }
00233         derived().finalize();
00234       }
00235       else
00236       {
00237         assignGeneric(other);
00238       }
00239       return derived();
00240     }
00241 
00242     template<typename OtherDerived>
00243     inline void assignGeneric(const OtherDerived& other)
00244     {
00245       //const bool transpose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit);
00246       eigen_assert(( ((internal::traits<Derived>::SupportedAccessPatterns&OuterRandomAccessPattern)==OuterRandomAccessPattern) ||
00247                   (!((Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit)))) &&
00248                   "the transpose operation is supposed to be handled in SparseMatrix::operator=");
00249 
00250       enum { Flip = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit) };
00251 
00252       const Index outerSize = other.outerSize();
00253       //typedef typename internal::conditional<transpose, LinkedVectorMatrix<Scalar,Flags&RowMajorBit>, Derived>::type TempType;
00254       // thanks to shallow copies, we always eval to a tempary
00255       Derived temp(other.rows(), other.cols());
00256 
00257       temp.reserve((std::max)(this->rows(),this->cols())*2);
00258       for (Index j=0; j<outerSize; ++j)
00259       {
00260         temp.startVec(j);
00261         for (typename OtherDerived::InnerIterator it(other.derived(), j); it; ++it)
00262         {
00263           Scalar v = it.value();
00264           temp.insertBackByOuterInner(Flip?it.index():j,Flip?j:it.index()) = v;
00265         }
00266       }
00267       temp.finalize();
00268 
00269       derived() = temp.markAsRValue();
00270     }
00271 
00272   public:
00273 
00274     template<typename Lhs, typename Rhs>
00275     inline Derived& operator=(const SparseSparseProduct<Lhs,Rhs>& product);
00276 
00277     friend std::ostream & operator << (std::ostream & s, const SparseMatrixBase& m)
00278     {
00279       typedef typename Derived::Nested Nested;
00280       typedef typename internal::remove_all<Nested>::type NestedCleaned;
00281 
00282       if (Flags&RowMajorBit)
00283       {
00284         const Nested nm(m.derived());
00285         for (Index row=0; row<nm.outerSize(); ++row)
00286         {
00287           Index col = 0;
00288           for (typename NestedCleaned::InnerIterator it(nm.derived(), row); it; ++it)
00289           {
00290             for ( ; col<it.index(); ++col)
00291               s << "0 ";
00292             s << it.value() << " ";
00293             ++col;
00294           }
00295           for ( ; col<m.cols(); ++col)
00296             s << "0 ";
00297           s << std::endl;
00298         }
00299       }
00300       else
00301       {
00302         const Nested nm(m.derived());
00303         if (m.cols() == 1) {
00304           Index row = 0;
00305           for (typename NestedCleaned::InnerIterator it(nm.derived(), 0); it; ++it)
00306           {
00307             for ( ; row<it.index(); ++row)
00308               s << "0" << std::endl;
00309             s << it.value() << std::endl;
00310             ++row;
00311           }
00312           for ( ; row<m.rows(); ++row)
00313             s << "0" << std::endl;
00314         }
00315         else
00316         {
00317           SparseMatrix<Scalar, RowMajorBit> trans = m;
00318           s << static_cast<const SparseMatrixBase<SparseMatrix<Scalar, RowMajorBit> >&>(trans);
00319         }
00320       }
00321       return s;
00322     }
00323 
00324     template<typename OtherDerived>
00325     Derived& operator+=(const SparseMatrixBase<OtherDerived>& other);
00326     template<typename OtherDerived>
00327     Derived& operator-=(const SparseMatrixBase<OtherDerived>& other);
00328 
00329     Derived& operator*=(const Scalar& other);
00330     Derived& operator/=(const Scalar& other);
00331 
00332     #define EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE \
00333       CwiseBinaryOp< \
00334         internal::scalar_product_op< \
00335           typename internal::scalar_product_traits< \
00336             typename internal::traits<Derived>::Scalar, \
00337             typename internal::traits<OtherDerived>::Scalar \
00338           >::ReturnType \
00339         >, \
00340         Derived, \
00341         OtherDerived \
00342       >
00343 
00344     template<typename OtherDerived>
00345     EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE
00346     cwiseProduct(const MatrixBase<OtherDerived> &other) const;
00347 
00348     // sparse * sparse
00349     template<typename OtherDerived>
00350     const typename SparseSparseProductReturnType<Derived,OtherDerived>::Type
00351     operator*(const SparseMatrixBase<OtherDerived> &other) const;
00352 
00353     // sparse * diagonal
00354     template<typename OtherDerived>
00355     const SparseDiagonalProduct<Derived,OtherDerived>
00356     operator*(const DiagonalBase<OtherDerived> &other) const;
00357 
00358     // diagonal * sparse
00359     template<typename OtherDerived> friend
00360     const SparseDiagonalProduct<OtherDerived,Derived>
00361     operator*(const DiagonalBase<OtherDerived> &lhs, const SparseMatrixBase& rhs)
00362     { return SparseDiagonalProduct<OtherDerived,Derived>(lhs.derived(), rhs.derived()); }
00363 
00365     template<typename OtherDerived> friend
00366     const typename DenseSparseProductReturnType<OtherDerived,Derived>::Type
00367     operator*(const MatrixBase<OtherDerived>& lhs, const Derived& rhs)
00368     { return typename DenseSparseProductReturnType<OtherDerived,Derived>::Type(lhs.derived(),rhs); }
00369 
00371     template<typename OtherDerived>
00372     const typename SparseDenseProductReturnType<Derived,OtherDerived>::Type
00373     operator*(const MatrixBase<OtherDerived> &other) const;
00374     
00376     SparseSymmetricPermutationProduct<Derived,Upper|Lower> twistedBy(const PermutationMatrix<Dynamic,Dynamic,Index>& perm) const
00377     {
00378       return SparseSymmetricPermutationProduct<Derived,Upper|Lower>(derived(), perm);
00379     }
00380 
00381     template<typename OtherDerived>
00382     Derived& operator*=(const SparseMatrixBase<OtherDerived>& other);
00383 
00384     #ifdef EIGEN2_SUPPORT
00385     // deprecated
00386     template<typename OtherDerived>
00387     typename internal::plain_matrix_type_column_major<OtherDerived>::type
00388     solveTriangular(const MatrixBase<OtherDerived>& other) const;
00389 
00390     // deprecated
00391     template<typename OtherDerived>
00392     void solveTriangularInPlace(MatrixBase<OtherDerived>& other) const;
00393     #endif // EIGEN2_SUPPORT
00394 
00395     template<int Mode>
00396     inline const SparseTriangularView<Derived, Mode> triangularView() const;
00397 
00398     template<unsigned int UpLo> inline const SparseSelfAdjointView<Derived, UpLo> selfadjointView() const;
00399     template<unsigned int UpLo> inline SparseSelfAdjointView<Derived, UpLo> selfadjointView();
00400 
00401     template<typename OtherDerived> Scalar dot(const MatrixBase<OtherDerived>& other) const;
00402     template<typename OtherDerived> Scalar dot(const SparseMatrixBase<OtherDerived>& other) const;
00403     RealScalar squaredNorm() const;
00404     RealScalar norm()  const;
00405 
00406     Transpose<Derived> transpose() { return derived(); }
00407     const Transpose<const Derived> transpose() const { return derived(); }
00408     const AdjointReturnType adjoint() const { return transpose(); }
00409 
00410     // sub-vector
00411     SparseInnerVectorSet<Derived,1> row(Index i);
00412     const SparseInnerVectorSet<Derived,1> row(Index i) const;
00413     SparseInnerVectorSet<Derived,1> col(Index j);
00414     const SparseInnerVectorSet<Derived,1> col(Index j) const;
00415     SparseInnerVectorSet<Derived,1> innerVector(Index outer);
00416     const SparseInnerVectorSet<Derived,1> innerVector(Index outer) const;
00417 
00418     // set of sub-vectors
00419     SparseInnerVectorSet<Derived,Dynamic> subrows(Index start, Index size);
00420     const SparseInnerVectorSet<Derived,Dynamic> subrows(Index start, Index size) const;
00421     SparseInnerVectorSet<Derived,Dynamic> subcols(Index start, Index size);
00422     const SparseInnerVectorSet<Derived,Dynamic> subcols(Index start, Index size) const;
00423     
00424     SparseInnerVectorSet<Derived,Dynamic> middleRows(Index start, Index size);
00425     const SparseInnerVectorSet<Derived,Dynamic> middleRows(Index start, Index size) const;
00426     SparseInnerVectorSet<Derived,Dynamic> middleCols(Index start, Index size);
00427     const SparseInnerVectorSet<Derived,Dynamic> middleCols(Index start, Index size) const;
00428     SparseInnerVectorSet<Derived,Dynamic> innerVectors(Index outerStart, Index outerSize);
00429     const SparseInnerVectorSet<Derived,Dynamic> innerVectors(Index outerStart, Index outerSize) const;
00430 
00432       template<typename DenseDerived>
00433       void evalTo(MatrixBase<DenseDerived>& dst) const
00434       {
00435         dst.setZero();
00436         for (Index j=0; j<outerSize(); ++j)
00437           for (typename Derived::InnerIterator i(derived(),j); i; ++i)
00438             dst.coeffRef(i.row(),i.col()) = i.value();
00439       }
00440 
00441       Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime> toDense() const
00442       {
00443         return derived();
00444       }
00445 
00446     template<typename OtherDerived>
00447     bool isApprox(const SparseMatrixBase<OtherDerived>& other,
00448                   RealScalar prec = NumTraits<Scalar>::dummy_precision()) const
00449     { return toDense().isApprox(other.toDense(),prec); }
00450 
00451     template<typename OtherDerived>
00452     bool isApprox(const MatrixBase<OtherDerived>& other,
00453                   RealScalar prec = NumTraits<Scalar>::dummy_precision()) const
00454     { return toDense().isApprox(other,prec); }
00455 
00461     inline const typename internal::eval<Derived>::type eval() const
00462     { return typename internal::eval<Derived>::type(derived()); }
00463 
00464     Scalar sum() const;
00465 
00466   protected:
00467 
00468     bool m_isRValue;
00469 };
00470 
00471 } // end namespace Eigen
00472 
00473 #endif // EIGEN_SPARSEMATRIXBASE_H