SparseProduct.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-2010 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_SPARSEPRODUCT_H
00026 #define EIGEN_SPARSEPRODUCT_H
00027 
00028 namespace Eigen { 
00029 
00030 template<typename Lhs, typename Rhs>
00031 struct SparseSparseProductReturnType
00032 {
00033   typedef typename internal::traits<Lhs>::Scalar Scalar;
00034   enum {
00035     LhsRowMajor = internal::traits<Lhs>::Flags & RowMajorBit,
00036     RhsRowMajor = internal::traits<Rhs>::Flags & RowMajorBit,
00037     TransposeRhs = (!LhsRowMajor) && RhsRowMajor,
00038     TransposeLhs = LhsRowMajor && (!RhsRowMajor)
00039   };
00040 
00041   typedef typename internal::conditional<TransposeLhs,
00042     SparseMatrix<Scalar,0>,
00043     typename internal::nested<Lhs,Rhs::RowsAtCompileTime>::type>::type LhsNested;
00044 
00045   typedef typename internal::conditional<TransposeRhs,
00046     SparseMatrix<Scalar,0>,
00047     typename internal::nested<Rhs,Lhs::RowsAtCompileTime>::type>::type RhsNested;
00048 
00049   typedef SparseSparseProduct<LhsNested, RhsNested> Type;
00050 };
00051 
00052 namespace internal {
00053 template<typename LhsNested, typename RhsNested>
00054 struct traits<SparseSparseProduct<LhsNested, RhsNested> >
00055 {
00056   typedef MatrixXpr XprKind;
00057   // clean the nested types:
00058   typedef typename remove_all<LhsNested>::type _LhsNested;
00059   typedef typename remove_all<RhsNested>::type _RhsNested;
00060   typedef typename _LhsNested::Scalar Scalar;
00061   typedef typename promote_index_type<typename traits<_LhsNested>::Index,
00062                                          typename traits<_RhsNested>::Index>::type Index;
00063 
00064   enum {
00065     LhsCoeffReadCost = _LhsNested::CoeffReadCost,
00066     RhsCoeffReadCost = _RhsNested::CoeffReadCost,
00067     LhsFlags = _LhsNested::Flags,
00068     RhsFlags = _RhsNested::Flags,
00069 
00070     RowsAtCompileTime    = _LhsNested::RowsAtCompileTime,
00071     ColsAtCompileTime    = _RhsNested::ColsAtCompileTime,
00072     MaxRowsAtCompileTime = _LhsNested::MaxRowsAtCompileTime,
00073     MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime,
00074 
00075     InnerSize = EIGEN_SIZE_MIN_PREFER_FIXED(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime),
00076 
00077     EvalToRowMajor = (RhsFlags & LhsFlags & RowMajorBit),
00078 
00079     RemovedBits = ~(EvalToRowMajor ? 0 : RowMajorBit),
00080 
00081     Flags = (int(LhsFlags | RhsFlags) & HereditaryBits & RemovedBits)
00082           | EvalBeforeAssigningBit
00083           | EvalBeforeNestingBit,
00084 
00085     CoeffReadCost = Dynamic
00086   };
00087 
00088   typedef Sparse StorageKind;
00089 };
00090 
00091 } // end namespace internal
00092 
00093 template<typename LhsNested, typename RhsNested>
00094 class SparseSparseProduct : internal::no_assignment_operator,
00095   public SparseMatrixBase<SparseSparseProduct<LhsNested, RhsNested> >
00096 {
00097   public:
00098 
00099     typedef SparseMatrixBase<SparseSparseProduct> Base;
00100     EIGEN_DENSE_PUBLIC_INTERFACE(SparseSparseProduct)
00101 
00102   private:
00103 
00104     typedef typename internal::traits<SparseSparseProduct>::_LhsNested _LhsNested;
00105     typedef typename internal::traits<SparseSparseProduct>::_RhsNested _RhsNested;
00106 
00107   public:
00108 
00109     template<typename Lhs, typename Rhs>
00110     EIGEN_STRONG_INLINE SparseSparseProduct(const Lhs& lhs, const Rhs& rhs)
00111       : m_lhs(lhs), m_rhs(rhs), m_tolerance(0), m_conservative(true)
00112     {
00113       init();
00114     }
00115 
00116     template<typename Lhs, typename Rhs>
00117     EIGEN_STRONG_INLINE SparseSparseProduct(const Lhs& lhs, const Rhs& rhs, RealScalar tolerance)
00118       : m_lhs(lhs), m_rhs(rhs), m_tolerance(tolerance), m_conservative(false)
00119     {
00120       init();
00121     }
00122 
00123     SparseSparseProduct pruned(Scalar reference = 0, RealScalar epsilon = NumTraits<RealScalar>::dummy_precision()) const
00124     {
00125       return SparseSparseProduct(m_lhs,m_rhs,internal::abs(reference)*epsilon);
00126     }
00127 
00128     template<typename Dest>
00129     void evalTo(Dest& result) const
00130     {
00131       if(m_conservative)
00132         internal::conservative_sparse_sparse_product_selector<_LhsNested, _RhsNested, Dest>::run(lhs(),rhs(),result);
00133       else
00134         internal::sparse_sparse_product_with_pruning_selector<_LhsNested, _RhsNested, Dest>::run(lhs(),rhs(),result,m_tolerance);
00135     }
00136 
00137     EIGEN_STRONG_INLINE Index rows() const { return m_lhs.rows(); }
00138     EIGEN_STRONG_INLINE Index cols() const { return m_rhs.cols(); }
00139 
00140     EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; }
00141     EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; }
00142 
00143   protected:
00144     void init()
00145     {
00146       eigen_assert(m_lhs.cols() == m_rhs.rows());
00147 
00148       enum {
00149         ProductIsValid = _LhsNested::ColsAtCompileTime==Dynamic
00150                       || _RhsNested::RowsAtCompileTime==Dynamic
00151                       || int(_LhsNested::ColsAtCompileTime)==int(_RhsNested::RowsAtCompileTime),
00152         AreVectors = _LhsNested::IsVectorAtCompileTime && _RhsNested::IsVectorAtCompileTime,
00153         SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(_LhsNested,_RhsNested)
00154       };
00155       // note to the lost user:
00156       //    * for a dot product use: v1.dot(v2)
00157       //    * for a coeff-wise product use: v1.cwise()*v2
00158       EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes),
00159         INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS)
00160       EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors),
00161         INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION)
00162       EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)
00163     }
00164 
00165     LhsNested m_lhs;
00166     RhsNested m_rhs;
00167     RealScalar m_tolerance;
00168     bool m_conservative;
00169 };
00170 
00171 // sparse = sparse * sparse
00172 template<typename Derived>
00173 template<typename Lhs, typename Rhs>
00174 inline Derived& SparseMatrixBase<Derived>::operator=(const SparseSparseProduct<Lhs,Rhs>& product)
00175 {
00176   product.evalTo(derived());
00177   return derived();
00178 }
00179 
00191 template<typename Derived>
00192 template<typename OtherDerived>
00193 inline const typename SparseSparseProductReturnType<Derived,OtherDerived>::Type
00194 SparseMatrixBase<Derived>::operator*(const SparseMatrixBase<OtherDerived> &other) const
00195 {
00196   return typename SparseSparseProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
00197 }
00198 
00199 } // end namespace Eigen
00200 
00201 #endif // EIGEN_SPARSEPRODUCT_H