1aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt// This file is part of Eigen, a lightweight C++ template library 2aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt// for linear algebra. 3aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt// 4aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> 5aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt// 6aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt// This Source Code Form is subject to the terms of the Mozilla 7aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt// Public License v. 2.0. If a copy of the MPL was not distributed 8aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt 102e0ec4fa0fb5162c441cd666f55fe76777e40d5eJan Engelhardt#ifndef EIGEN_SPARSE_DIAGONAL_PRODUCT_H 11aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt#define EIGEN_SPARSE_DIAGONAL_PRODUCT_H 12aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt 13aa37acc1423126f555135935c687eb91995b9440Jan Engelhardtnamespace Eigen { 14aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt 15aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt// The product of a diagonal matrix with a sparse matrix can be easily 16aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt// implemented using expression template. 17aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt// We have two consider very different cases: 18aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt// 1 - diag * row-major sparse 19aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt// => each inner vector <=> scalar * sparse vector product 2041a4cea0f4109fb76762dca073c3c1217658ee06Jan Engelhardt// => so we can reuse CwiseUnaryOp::InnerIterator 21aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt// 2 - diag * col-major sparse 2261cc52b6f9edfa3efb1d0c9ea9531abb42828ec2Jan Engelhardt// => each inner vector <=> densevector * sparse vector cwise product 23aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt// => again, we can reuse specialization of CwiseBinaryOp::InnerIterator 24aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt// for that particular case 2561cc52b6f9edfa3efb1d0c9ea9531abb42828ec2Jan Engelhardt// The two other cases are symmetric. 2661cc52b6f9edfa3efb1d0c9ea9531abb42828ec2Jan Engelhardt 2761cc52b6f9edfa3efb1d0c9ea9531abb42828ec2Jan Engelhardtnamespace internal { 28aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt 29aa37acc1423126f555135935c687eb91995b9440Jan Engelhardttemplate<typename Lhs, typename Rhs> 30aa37acc1423126f555135935c687eb91995b9440Jan Engelhardtstruct traits<SparseDiagonalProduct<Lhs, Rhs> > 31aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt{ 32aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt typedef typename remove_all<Lhs>::type _Lhs; 3341a4cea0f4109fb76762dca073c3c1217658ee06Jan Engelhardt typedef typename remove_all<Rhs>::type _Rhs; 3441a4cea0f4109fb76762dca073c3c1217658ee06Jan Engelhardt typedef typename _Lhs::Scalar Scalar; 3541a4cea0f4109fb76762dca073c3c1217658ee06Jan Engelhardt typedef typename promote_index_type<typename traits<Lhs>::Index, 3641a4cea0f4109fb76762dca073c3c1217658ee06Jan Engelhardt typename traits<Rhs>::Index>::type Index; 3741a4cea0f4109fb76762dca073c3c1217658ee06Jan Engelhardt typedef Sparse StorageKind; 3841a4cea0f4109fb76762dca073c3c1217658ee06Jan Engelhardt typedef MatrixXpr XprKind; 3941a4cea0f4109fb76762dca073c3c1217658ee06Jan Engelhardt enum { 4061cc52b6f9edfa3efb1d0c9ea9531abb42828ec2Jan Engelhardt RowsAtCompileTime = _Lhs::RowsAtCompileTime, 4161cc52b6f9edfa3efb1d0c9ea9531abb42828ec2Jan Engelhardt ColsAtCompileTime = _Rhs::ColsAtCompileTime, 4261cc52b6f9edfa3efb1d0c9ea9531abb42828ec2Jan Engelhardt 4361cc52b6f9edfa3efb1d0c9ea9531abb42828ec2Jan Engelhardt MaxRowsAtCompileTime = _Lhs::MaxRowsAtCompileTime, 4441a4cea0f4109fb76762dca073c3c1217658ee06Jan Engelhardt MaxColsAtCompileTime = _Rhs::MaxColsAtCompileTime, 45aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt 46aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt SparseFlags = is_diagonal<_Lhs>::ret ? int(_Rhs::Flags) : int(_Lhs::Flags), 47aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt Flags = (SparseFlags&RowMajorBit), 48aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt CoeffReadCost = Dynamic 49aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt }; 50aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt}; 51aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt 52aa37acc1423126f555135935c687eb91995b9440Jan Engelhardtenum {SDP_IsDiagonal, SDP_IsSparseRowMajor, SDP_IsSparseColMajor}; 53aa37acc1423126f555135935c687eb91995b9440Jan Engelhardttemplate<typename Lhs, typename Rhs, typename SparseDiagonalProductType, int RhsMode, int LhsMode> 54aa37acc1423126f555135935c687eb91995b9440Jan Engelhardtclass sparse_diagonal_product_inner_iterator_selector; 55aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt 56aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt} // end namespace internal 57aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt 58aa37acc1423126f555135935c687eb91995b9440Jan Engelhardttemplate<typename Lhs, typename Rhs> 59aa37acc1423126f555135935c687eb91995b9440Jan Engelhardtclass SparseDiagonalProduct 60aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt : public SparseMatrixBase<SparseDiagonalProduct<Lhs,Rhs> >, 61aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt internal::no_assignment_operator 62aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt{ 63aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt typedef typename Lhs::Nested LhsNested; 64aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt typedef typename Rhs::Nested RhsNested; 65aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt 66aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt typedef typename internal::remove_all<LhsNested>::type _LhsNested; 67aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt typedef typename internal::remove_all<RhsNested>::type _RhsNested; 68aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt 69aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt enum { 70aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt LhsMode = internal::is_diagonal<_LhsNested>::ret ? internal::SDP_IsDiagonal 71aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt : (_LhsNested::Flags&RowMajorBit) ? internal::SDP_IsSparseRowMajor : internal::SDP_IsSparseColMajor, 72aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt RhsMode = internal::is_diagonal<_RhsNested>::ret ? internal::SDP_IsDiagonal 73aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt : (_RhsNested::Flags&RowMajorBit) ? internal::SDP_IsSparseRowMajor : internal::SDP_IsSparseColMajor 74aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt }; 75aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt 76aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt public: 77aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt 78aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt EIGEN_SPARSE_PUBLIC_INTERFACE(SparseDiagonalProduct) 79aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt 80aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt typedef internal::sparse_diagonal_product_inner_iterator_selector 81aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt <_LhsNested,_RhsNested,SparseDiagonalProduct,LhsMode,RhsMode> InnerIterator; 82aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt 83aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt // We do not want ReverseInnerIterator for diagonal-sparse products, 84aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt // but this dummy declaration is needed to make diag * sparse * diag compile. 85aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt class ReverseInnerIterator; 86aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt 87aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt EIGEN_STRONG_INLINE SparseDiagonalProduct(const Lhs& lhs, const Rhs& rhs) 88aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt : m_lhs(lhs), m_rhs(rhs) 89aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt { 90aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt eigen_assert(lhs.cols() == rhs.rows() && "invalid sparse matrix * diagonal matrix product"); 91aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt } 92aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt 93aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt EIGEN_STRONG_INLINE Index rows() const { return m_lhs.rows(); } 94aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt EIGEN_STRONG_INLINE Index cols() const { return m_rhs.cols(); } 95aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt 96aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; } 97aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; } 98aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt 99aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt protected: 100aa37acc1423126f555135935c687eb91995b9440Jan Engelhardt LhsNested m_lhs; 101a93142d5f55db74ebd7d49be9bd88f7a499ded40Jan Engelhardt RhsNested m_rhs; 102a93142d5f55db74ebd7d49be9bd88f7a499ded40Jan Engelhardt}; 103a93142d5f55db74ebd7d49be9bd88f7a499ded40Jan Engelhardt 104a93142d5f55db74ebd7d49be9bd88f7a499ded40Jan Engelhardtnamespace internal { 105a93142d5f55db74ebd7d49be9bd88f7a499ded40Jan Engelhardt 106a93142d5f55db74ebd7d49be9bd88f7a499ded40Jan Engelhardttemplate<typename Lhs, typename Rhs, typename SparseDiagonalProductType> 1078b5bdea659f1fb86b3288a2568ab104a90b914e5Jan Engelhardtclass sparse_diagonal_product_inner_iterator_selector 108a93142d5f55db74ebd7d49be9bd88f7a499ded40Jan Engelhardt<Lhs,Rhs,SparseDiagonalProductType,SDP_IsDiagonal,SDP_IsSparseRowMajor> 109a93142d5f55db74ebd7d49be9bd88f7a499ded40Jan Engelhardt : public CwiseUnaryOp<scalar_multiple_op<typename Lhs::Scalar>,const Rhs>::InnerIterator 110dfe99f1bf291b4b954d3608dbe95a43e16a8bb49Jan Engelhardt{ 111dfe99f1bf291b4b954d3608dbe95a43e16a8bb49Jan Engelhardt typedef typename CwiseUnaryOp<scalar_multiple_op<typename Lhs::Scalar>,const Rhs>::InnerIterator Base; 1120eff54bd407aae6b99c3b189d356929e399b5a38Jan Engelhardt typedef typename Lhs::Index Index; 1130eff54bd407aae6b99c3b189d356929e399b5a38Jan Engelhardt public: 1148b5bdea659f1fb86b3288a2568ab104a90b914e5Jan Engelhardt inline sparse_diagonal_product_inner_iterator_selector( 1158b5bdea659f1fb86b3288a2568ab104a90b914e5Jan Engelhardt const SparseDiagonalProductType& expr, Index outer) 116d78254d7f9d18ef76377a3013302430cce8ea702Jan Engelhardt : Base(expr.rhs()*(expr.lhs().diagonal().coeff(outer)), outer) 117d78254d7f9d18ef76377a3013302430cce8ea702Jan Engelhardt {} 118d78254d7f9d18ef76377a3013302430cce8ea702Jan Engelhardt}; 119d78254d7f9d18ef76377a3013302430cce8ea702Jan Engelhardt 120d78254d7f9d18ef76377a3013302430cce8ea702Jan Engelhardttemplate<typename Lhs, typename Rhs, typename SparseDiagonalProductType> 121a93142d5f55db74ebd7d49be9bd88f7a499ded40Jan Engelhardtclass sparse_diagonal_product_inner_iterator_selector 122a93142d5f55db74ebd7d49be9bd88f7a499ded40Jan Engelhardt<Lhs,Rhs,SparseDiagonalProductType,SDP_IsDiagonal,SDP_IsSparseColMajor> 123a93142d5f55db74ebd7d49be9bd88f7a499ded40Jan Engelhardt : public CwiseBinaryOp< 1248b5bdea659f1fb86b3288a2568ab104a90b914e5Jan Engelhardt scalar_product_op<typename Lhs::Scalar>, 125a93142d5f55db74ebd7d49be9bd88f7a499ded40Jan Engelhardt const typename Rhs::ConstInnerVectorReturnType, 126a93142d5f55db74ebd7d49be9bd88f7a499ded40Jan Engelhardt const typename Lhs::DiagonalVectorType>::InnerIterator 127dfe99f1bf291b4b954d3608dbe95a43e16a8bb49Jan Engelhardt{ 128dfe99f1bf291b4b954d3608dbe95a43e16a8bb49Jan Engelhardt typedef typename CwiseBinaryOp< 129dfe99f1bf291b4b954d3608dbe95a43e16a8bb49Jan Engelhardt scalar_product_op<typename Lhs::Scalar>, 130dfe99f1bf291b4b954d3608dbe95a43e16a8bb49Jan Engelhardt const typename Rhs::ConstInnerVectorReturnType, 1310eff54bd407aae6b99c3b189d356929e399b5a38Jan Engelhardt const typename Lhs::DiagonalVectorType>::InnerIterator Base; 1320eff54bd407aae6b99c3b189d356929e399b5a38Jan Engelhardt typedef typename Lhs::Index Index; 1330eff54bd407aae6b99c3b189d356929e399b5a38Jan Engelhardt Index m_outer; 1340eff54bd407aae6b99c3b189d356929e399b5a38Jan Engelhardt public: 135dfe99f1bf291b4b954d3608dbe95a43e16a8bb49Jan Engelhardt inline sparse_diagonal_product_inner_iterator_selector( 136a93142d5f55db74ebd7d49be9bd88f7a499ded40Jan Engelhardt const SparseDiagonalProductType& expr, Index outer) 137a93142d5f55db74ebd7d49be9bd88f7a499ded40Jan Engelhardt : Base(expr.rhs().innerVector(outer) .cwiseProduct(expr.lhs().diagonal()), 0), m_outer(outer) 138a93142d5f55db74ebd7d49be9bd88f7a499ded40Jan Engelhardt {} 1398b5bdea659f1fb86b3288a2568ab104a90b914e5Jan Engelhardt 1408b5bdea659f1fb86b3288a2568ab104a90b914e5Jan Engelhardt inline Index outer() const { return m_outer; } 1418b5bdea659f1fb86b3288a2568ab104a90b914e5Jan Engelhardt inline Index col() const { return m_outer; } 1428b5bdea659f1fb86b3288a2568ab104a90b914e5Jan Engelhardt}; 143a93142d5f55db74ebd7d49be9bd88f7a499ded40Jan Engelhardt 144a93142d5f55db74ebd7d49be9bd88f7a499ded40Jan Engelhardttemplate<typename Lhs, typename Rhs, typename SparseDiagonalProductType> 145a93142d5f55db74ebd7d49be9bd88f7a499ded40Jan Engelhardtclass sparse_diagonal_product_inner_iterator_selector 14604bb988275ac76815a15788a7fc75ac78f3bb833Jan Engelhardt<Lhs,Rhs,SparseDiagonalProductType,SDP_IsSparseColMajor,SDP_IsDiagonal> 147f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardt : public CwiseUnaryOp<scalar_multiple_op<typename Rhs::Scalar>,const Lhs>::InnerIterator 148f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardt{ 149f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardt typedef typename CwiseUnaryOp<scalar_multiple_op<typename Rhs::Scalar>,const Lhs>::InnerIterator Base; 150f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardt typedef typename Lhs::Index Index; 151f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardt public: 152f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardt inline sparse_diagonal_product_inner_iterator_selector( 153f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardt const SparseDiagonalProductType& expr, Index outer) 154f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardt : Base(expr.lhs()*expr.rhs().diagonal().coeff(outer), outer) 155f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardt {} 156f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardt}; 157f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardt 158f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardttemplate<typename Lhs, typename Rhs, typename SparseDiagonalProductType> 159f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardtclass sparse_diagonal_product_inner_iterator_selector 160f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardt<Lhs,Rhs,SparseDiagonalProductType,SDP_IsSparseRowMajor,SDP_IsDiagonal> 161f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardt : public CwiseBinaryOp< 162f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardt scalar_product_op<typename Rhs::Scalar>, 163f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardt const typename Lhs::ConstInnerVectorReturnType, 164f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardt const Transpose<const typename Rhs::DiagonalVectorType> >::InnerIterator 165f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardt{ 166f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardt typedef typename CwiseBinaryOp< 167f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardt scalar_product_op<typename Rhs::Scalar>, 168f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardt const typename Lhs::ConstInnerVectorReturnType, 169f012b3c9190cd95ac170072f759a97575613ea07Jan Engelhardt const Transpose<const typename Rhs::DiagonalVectorType> >::InnerIterator Base; 17004bb988275ac76815a15788a7fc75ac78f3bb833Jan Engelhardt typedef typename Lhs::Index Index; 17104bb988275ac76815a15788a7fc75ac78f3bb833Jan Engelhardt Index m_outer; 17204bb988275ac76815a15788a7fc75ac78f3bb833Jan Engelhardt public: 17304bb988275ac76815a15788a7fc75ac78f3bb833Jan Engelhardt inline sparse_diagonal_product_inner_iterator_selector( 17404bb988275ac76815a15788a7fc75ac78f3bb833Jan Engelhardt const SparseDiagonalProductType& expr, Index outer) 17504bb988275ac76815a15788a7fc75ac78f3bb833Jan Engelhardt : Base(expr.lhs().innerVector(outer) .cwiseProduct(expr.rhs().diagonal().transpose()), 0), m_outer(outer) 17604bb988275ac76815a15788a7fc75ac78f3bb833Jan Engelhardt {} 17704bb988275ac76815a15788a7fc75ac78f3bb833Jan Engelhardt 17804bb988275ac76815a15788a7fc75ac78f3bb833Jan Engelhardt inline Index outer() const { return m_outer; } 17904bb988275ac76815a15788a7fc75ac78f3bb833Jan Engelhardt inline Index row() const { return m_outer; } 18004bb988275ac76815a15788a7fc75ac78f3bb833Jan Engelhardt}; 18104bb988275ac76815a15788a7fc75ac78f3bb833Jan Engelhardt 182564eaf48e14411803a353206eefbb89d525c63ffJan Engelhardt} // end namespace internal 183564eaf48e14411803a353206eefbb89d525c63ffJan Engelhardt 18404bb988275ac76815a15788a7fc75ac78f3bb833Jan Engelhardt// SparseMatrixBase functions 18504bb988275ac76815a15788a7fc75ac78f3bb833Jan Engelhardt 18604bb988275ac76815a15788a7fc75ac78f3bb833Jan Engelhardttemplate<typename Derived> 18704bb988275ac76815a15788a7fc75ac78f3bb833Jan Engelhardttemplate<typename OtherDerived> 1888bf513ada0aae0e4b1ac5160113fc532c2f525d0Jan Engelhardtconst SparseDiagonalProduct<Derived,OtherDerived> 1898bf513ada0aae0e4b1ac5160113fc532c2f525d0Jan EngelhardtSparseMatrixBase<Derived>::operator*(const DiagonalBase<OtherDerived> &other) const 1908bf513ada0aae0e4b1ac5160113fc532c2f525d0Jan Engelhardt{ 191564eaf48e14411803a353206eefbb89d525c63ffJan Engelhardt return SparseDiagonalProduct<Derived,OtherDerived>(this->derived(), other.derived()); 192bc438c4cbdab09fafbbceecddd54e44e4234a4a1Jan Engelhardt} 193bc438c4cbdab09fafbbceecddd54e44e4234a4a1Jan Engelhardt 194564eaf48e14411803a353206eefbb89d525c63ffJan Engelhardt} // end namespace Eigen 19504bb988275ac76815a15788a7fc75ac78f3bb833Jan Engelhardt 19604bb988275ac76815a15788a7fc75ac78f3bb833Jan Engelhardt#endif // EIGEN_SPARSE_DIAGONAL_PRODUCT_H 197564eaf48e14411803a353206eefbb89d525c63ffJan Engelhardt