src/SparseCore/SparseProduct.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2008-2015 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_SPARSEPRODUCT_H
0011 #define EIGEN_SPARSEPRODUCT_H
0012
0013 namespace Eigen {
0014
0015 /** \returns an expression of the product of two sparse matrices.
0016   * By default a conservative product preserving the symbolic non zeros is performed.
0017   * The automatic pruning of the small values can be achieved by calling the pruned() function
0018   * in which case a totally different product algorithm is employed:
0019   * \code
0020   * C = (A*B).pruned();             // suppress numerical zeros (exact)
0021   * C = (A*B).pruned(ref);
0022   * C = (A*B).pruned(ref,epsilon);
0023   * \endcode
0024   * where \c ref is a meaningful non zero reference value.
0025   * */
0026 template<typename Derived>
0027 template<typename OtherDerived>
0028 inline const Product<Derived,OtherDerived,AliasFreeProduct>
0029 SparseMatrixBase<Derived>::operator*(const SparseMatrixBase<OtherDerived> &other) const
0030 {
0031   return Product<Derived,OtherDerived,AliasFreeProduct>(derived(), other.derived());
0032 }
0033
0034 namespace internal {
0035
0036 // sparse * sparse
0037 template<typename Lhs, typename Rhs, int ProductType>
0038 struct generic_product_impl<Lhs, Rhs, SparseShape, SparseShape, ProductType>
0039 {
0040   template<typename Dest>
0041   static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
0042   {
0043     evalTo(dst, lhs, rhs, typename evaluator_traits<Dest>::Shape());
0044   }
0045
0046   // dense += sparse * sparse
0047   template<typename Dest,typename ActualLhs>
0048   static void addTo(Dest& dst, const ActualLhs& lhs, const Rhs& rhs, typename enable_if<is_same<typename evaluator_traits<Dest>::Shape,DenseShape>::value,int*>::type* = 0)
0049   {
0050     typedef typename nested_eval<ActualLhs,Dynamic>::type LhsNested;
0051     typedef typename nested_eval<Rhs,Dynamic>::type RhsNested;
0052     LhsNested lhsNested(lhs);
0053     RhsNested rhsNested(rhs);
0054     internal::sparse_sparse_to_dense_product_selector<typename remove_all<LhsNested>::type,
0055                                                       typename remove_all<RhsNested>::type, Dest>::run(lhsNested,rhsNested,dst);
0056   }
0057
0058   // dense -= sparse * sparse
0059   template<typename Dest>
0060   static void subTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, typename enable_if<is_same<typename evaluator_traits<Dest>::Shape,DenseShape>::value,int*>::type* = 0)
0061   {
0062     addTo(dst, -lhs, rhs);
0063   }
0064
0065 protected:
0066
0067   // sparse = sparse * sparse
0068   template<typename Dest>
0069   static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, SparseShape)
0070   {
0071     typedef typename nested_eval<Lhs,Dynamic>::type LhsNested;
0072     typedef typename nested_eval<Rhs,Dynamic>::type RhsNested;
0073     LhsNested lhsNested(lhs);
0074     RhsNested rhsNested(rhs);
0075     internal::conservative_sparse_sparse_product_selector<typename remove_all<LhsNested>::type,
0076                                                           typename remove_all<RhsNested>::type, Dest>::run(lhsNested,rhsNested,dst);
0077   }
0078
0079   // dense = sparse * sparse
0080   template<typename Dest>
0081   static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, DenseShape)
0082   {
0083     dst.setZero();
0084     addTo(dst, lhs, rhs);
0085   }
0086 };
0087
0088 // sparse * sparse-triangular
0089 template<typename Lhs, typename Rhs, int ProductType>
0090 struct generic_product_impl<Lhs, Rhs, SparseShape, SparseTriangularShape, ProductType>
0091  : public generic_product_impl<Lhs, Rhs, SparseShape, SparseShape, ProductType>
0092 {};
0093
0094 // sparse-triangular * sparse
0095 template<typename Lhs, typename Rhs, int ProductType>
0096 struct generic_product_impl<Lhs, Rhs, SparseTriangularShape, SparseShape, ProductType>
0097  : public generic_product_impl<Lhs, Rhs, SparseShape, SparseShape, ProductType>
0098 {};
0099
0100 // dense = sparse-product (can be sparse*sparse, sparse*perm, etc.)
0101 template< typename DstXprType, typename Lhs, typename Rhs>
0102 struct Assignment<DstXprType, Product<Lhs,Rhs,AliasFreeProduct>, internal::assign_op<typename DstXprType::Scalar,typename Product<Lhs,Rhs,AliasFreeProduct>::Scalar>, Sparse2Dense>
0103 {
0104   typedef Product<Lhs,Rhs,AliasFreeProduct> SrcXprType;
0105   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &)
0106   {
0107     Index dstRows = src.rows();
0108     Index dstCols = src.cols();
0109     if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
0110       dst.resize(dstRows, dstCols);
0111
0112     generic_product_impl<Lhs, Rhs>::evalTo(dst,src.lhs(),src.rhs());
0113   }
0114 };
0115
0116 // dense += sparse-product (can be sparse*sparse, sparse*perm, etc.)
0117 template< typename DstXprType, typename Lhs, typename Rhs>
0118 struct Assignment<DstXprType, Product<Lhs,Rhs,AliasFreeProduct>, internal::add_assign_op<typename DstXprType::Scalar,typename Product<Lhs,Rhs,AliasFreeProduct>::Scalar>, Sparse2Dense>
0119 {
0120   typedef Product<Lhs,Rhs,AliasFreeProduct> SrcXprType;
0121   static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &)
0122   {
0123     generic_product_impl<Lhs, Rhs>::addTo(dst,src.lhs(),src.rhs());
0124   }
0125 };
0126
0127 // dense -= sparse-product (can be sparse*sparse, sparse*perm, etc.)
0128 template< typename DstXprType, typename Lhs, typename Rhs>
0129 struct Assignment<DstXprType, Product<Lhs,Rhs,AliasFreeProduct>, internal::sub_assign_op<typename DstXprType::Scalar,typename Product<Lhs,Rhs,AliasFreeProduct>::Scalar>, Sparse2Dense>
0130 {
0131   typedef Product<Lhs,Rhs,AliasFreeProduct> SrcXprType;
0132   static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar> &)
0133   {
0134     generic_product_impl<Lhs, Rhs>::subTo(dst,src.lhs(),src.rhs());
0135   }
0136 };
0137
0138 template<typename Lhs, typename Rhs, int Options>
0139 struct unary_evaluator<SparseView<Product<Lhs, Rhs, Options> >, IteratorBased>
0140  : public evaluator<typename Product<Lhs, Rhs, DefaultProduct>::PlainObject>
0141 {
0142   typedef SparseView<Product<Lhs, Rhs, Options> > XprType;
0143   typedef typename XprType::PlainObject PlainObject;
0144   typedef evaluator<PlainObject> Base;
0145
0146   explicit unary_evaluator(const XprType& xpr)
0147     : m_result(xpr.rows(), xpr.cols())
0148   {
0149     using std::abs;
0150     ::new (static_cast<Base*>(this)) Base(m_result);
0151     typedef typename nested_eval<Lhs,Dynamic>::type LhsNested;
0152     typedef typename nested_eval<Rhs,Dynamic>::type RhsNested;
0153     LhsNested lhsNested(xpr.nestedExpression().lhs());
0154     RhsNested rhsNested(xpr.nestedExpression().rhs());
0155
0156     internal::sparse_sparse_product_with_pruning_selector<typename remove_all<LhsNested>::type,
0157                                                           typename remove_all<RhsNested>::type, PlainObject>::run(lhsNested,rhsNested,m_result,
0158                                                                                                                   abs(xpr.reference())*xpr.epsilon());
0159   }
0160
0161 protected:
0162   PlainObject m_result;
0163 };
0164
0165 } // end namespace internal
0166
0167 // sparse matrix = sparse-product (can be sparse*sparse, sparse*perm, etc.)
0168 template<typename Scalar, int _Options, typename _StorageIndex>
0169 template<typename Lhs, typename Rhs>
0170 SparseMatrix<Scalar,_Options,_StorageIndex>& SparseMatrix<Scalar,_Options,_StorageIndex>::operator=(const Product<Lhs,Rhs,AliasFreeProduct>& src)
0171 {
0172   // std::cout << "in Assignment : " << DstOptions << "\n";
0173   SparseMatrix dst(src.rows(),src.cols());
0174   internal::generic_product_impl<Lhs, Rhs>::evalTo(dst,src.lhs(),src.rhs());
0175   this->swap(dst);
0176   return *this;
0177 }
0178
0179 } // end namespace Eigen
0180
0181 #endif // EIGEN_SPARSEPRODUCT_H