$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_SPARSE_DIAGONAL_PRODUCT_H 00011 #define EIGEN_SPARSE_DIAGONAL_PRODUCT_H 00012 00013 namespace Eigen { 00014 00015 // The product of a diagonal matrix with a sparse matrix can be easily 00016 // implemented using expression template. 00017 // We have two consider very different cases: 00018 // 1 - diag * row-major sparse 00019 // => each inner vector <=> scalar * sparse vector product 00020 // => so we can reuse CwiseUnaryOp::InnerIterator 00021 // 2 - diag * col-major sparse 00022 // => each inner vector <=> densevector * sparse vector cwise product 00023 // => again, we can reuse specialization of CwiseBinaryOp::InnerIterator 00024 // for that particular case 00025 // The two other cases are symmetric. 00026 00027 namespace internal { 00028 00029 template<typename Lhs, typename Rhs> 00030 struct traits<SparseDiagonalProduct<Lhs, Rhs> > 00031 { 00032 typedef typename remove_all<Lhs>::type _Lhs; 00033 typedef typename remove_all<Rhs>::type _Rhs; 00034 typedef typename _Lhs::Scalar Scalar; 00035 typedef typename promote_index_type<typename traits<Lhs>::Index, 00036 typename traits<Rhs>::Index>::type Index; 00037 typedef Sparse StorageKind; 00038 typedef MatrixXpr XprKind; 00039 enum { 00040 RowsAtCompileTime = _Lhs::RowsAtCompileTime, 00041 ColsAtCompileTime = _Rhs::ColsAtCompileTime, 00042 00043 MaxRowsAtCompileTime = _Lhs::MaxRowsAtCompileTime, 00044 MaxColsAtCompileTime = _Rhs::MaxColsAtCompileTime, 00045 00046 SparseFlags = is_diagonal<_Lhs>::ret ? int(_Rhs::Flags) : int(_Lhs::Flags), 00047 Flags = (SparseFlags&RowMajorBit), 00048 CoeffReadCost = Dynamic 00049 }; 00050 }; 00051 00052 enum {SDP_IsDiagonal, SDP_IsSparseRowMajor, SDP_IsSparseColMajor}; 00053 template<typename Lhs, typename Rhs, typename SparseDiagonalProductType, int RhsMode, int LhsMode> 00054 class sparse_diagonal_product_inner_iterator_selector; 00055 00056 } // end namespace internal 00057 00058 template<typename Lhs, typename Rhs> 00059 class SparseDiagonalProduct 00060 : public SparseMatrixBase<SparseDiagonalProduct<Lhs,Rhs> >, 00061 internal::no_assignment_operator 00062 { 00063 typedef typename Lhs::Nested LhsNested; 00064 typedef typename Rhs::Nested RhsNested; 00065 00066 typedef typename internal::remove_all<LhsNested>::type _LhsNested; 00067 typedef typename internal::remove_all<RhsNested>::type _RhsNested; 00068 00069 enum { 00070 LhsMode = internal::is_diagonal<_LhsNested>::ret ? internal::SDP_IsDiagonal 00071 : (_LhsNested::Flags&RowMajorBit) ? internal::SDP_IsSparseRowMajor : internal::SDP_IsSparseColMajor, 00072 RhsMode = internal::is_diagonal<_RhsNested>::ret ? internal::SDP_IsDiagonal 00073 : (_RhsNested::Flags&RowMajorBit) ? internal::SDP_IsSparseRowMajor : internal::SDP_IsSparseColMajor 00074 }; 00075 00076 public: 00077 00078 EIGEN_SPARSE_PUBLIC_INTERFACE(SparseDiagonalProduct) 00079 00080 typedef internal::sparse_diagonal_product_inner_iterator_selector 00081 <_LhsNested,_RhsNested,SparseDiagonalProduct,LhsMode,RhsMode> InnerIterator; 00082 00083 // We do not want ReverseInnerIterator for diagonal-sparse products, 00084 // but this dummy declaration is needed to make diag * sparse * diag compile. 00085 class ReverseInnerIterator; 00086 00087 EIGEN_STRONG_INLINE SparseDiagonalProduct(const Lhs& lhs, const Rhs& rhs) 00088 : m_lhs(lhs), m_rhs(rhs) 00089 { 00090 eigen_assert(lhs.cols() == rhs.rows() && "invalid sparse matrix * diagonal matrix product"); 00091 } 00092 00093 EIGEN_STRONG_INLINE Index rows() const { return m_lhs.rows(); } 00094 EIGEN_STRONG_INLINE Index cols() const { return m_rhs.cols(); } 00095 00096 EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; } 00097 EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; } 00098 00099 protected: 00100 LhsNested m_lhs; 00101 RhsNested m_rhs; 00102 }; 00103 00104 namespace internal { 00105 00106 template<typename Lhs, typename Rhs, typename SparseDiagonalProductType> 00107 class sparse_diagonal_product_inner_iterator_selector 00108 <Lhs,Rhs,SparseDiagonalProductType,SDP_IsDiagonal,SDP_IsSparseRowMajor> 00109 : public CwiseUnaryOp<scalar_multiple_op<typename Lhs::Scalar>,const Rhs>::InnerIterator 00110 { 00111 typedef typename CwiseUnaryOp<scalar_multiple_op<typename Lhs::Scalar>,const Rhs>::InnerIterator Base; 00112 typedef typename Lhs::Index Index; 00113 public: 00114 inline sparse_diagonal_product_inner_iterator_selector( 00115 const SparseDiagonalProductType& expr, Index outer) 00116 : Base(expr.rhs()*(expr.lhs().diagonal().coeff(outer)), outer) 00117 {} 00118 }; 00119 00120 template<typename Lhs, typename Rhs, typename SparseDiagonalProductType> 00121 class sparse_diagonal_product_inner_iterator_selector 00122 <Lhs,Rhs,SparseDiagonalProductType,SDP_IsDiagonal,SDP_IsSparseColMajor> 00123 : public CwiseBinaryOp< 00124 scalar_product_op<typename Lhs::Scalar>, 00125 const typename Rhs::ConstInnerVectorReturnType, 00126 const typename Lhs::DiagonalVectorType>::InnerIterator 00127 { 00128 typedef typename CwiseBinaryOp< 00129 scalar_product_op<typename Lhs::Scalar>, 00130 const typename Rhs::ConstInnerVectorReturnType, 00131 const typename Lhs::DiagonalVectorType>::InnerIterator Base; 00132 typedef typename Lhs::Index Index; 00133 Index m_outer; 00134 public: 00135 inline sparse_diagonal_product_inner_iterator_selector( 00136 const SparseDiagonalProductType& expr, Index outer) 00137 : Base(expr.rhs().innerVector(outer) .cwiseProduct(expr.lhs().diagonal()), 0), m_outer(outer) 00138 {} 00139 00140 inline Index outer() const { return m_outer; } 00141 inline Index col() const { return m_outer; } 00142 }; 00143 00144 template<typename Lhs, typename Rhs, typename SparseDiagonalProductType> 00145 class sparse_diagonal_product_inner_iterator_selector 00146 <Lhs,Rhs,SparseDiagonalProductType,SDP_IsSparseColMajor,SDP_IsDiagonal> 00147 : public CwiseUnaryOp<scalar_multiple_op<typename Rhs::Scalar>,const Lhs>::InnerIterator 00148 { 00149 typedef typename CwiseUnaryOp<scalar_multiple_op<typename Rhs::Scalar>,const Lhs>::InnerIterator Base; 00150 typedef typename Lhs::Index Index; 00151 public: 00152 inline sparse_diagonal_product_inner_iterator_selector( 00153 const SparseDiagonalProductType& expr, Index outer) 00154 : Base(expr.lhs()*expr.rhs().diagonal().coeff(outer), outer) 00155 {} 00156 }; 00157 00158 template<typename Lhs, typename Rhs, typename SparseDiagonalProductType> 00159 class sparse_diagonal_product_inner_iterator_selector 00160 <Lhs,Rhs,SparseDiagonalProductType,SDP_IsSparseRowMajor,SDP_IsDiagonal> 00161 : public CwiseBinaryOp< 00162 scalar_product_op<typename Rhs::Scalar>, 00163 const typename Lhs::ConstInnerVectorReturnType, 00164 const Transpose<const typename Rhs::DiagonalVectorType> >::InnerIterator 00165 { 00166 typedef typename CwiseBinaryOp< 00167 scalar_product_op<typename Rhs::Scalar>, 00168 const typename Lhs::ConstInnerVectorReturnType, 00169 const Transpose<const typename Rhs::DiagonalVectorType> >::InnerIterator Base; 00170 typedef typename Lhs::Index Index; 00171 Index m_outer; 00172 public: 00173 inline sparse_diagonal_product_inner_iterator_selector( 00174 const SparseDiagonalProductType& expr, Index outer) 00175 : Base(expr.lhs().innerVector(outer) .cwiseProduct(expr.rhs().diagonal().transpose()), 0), m_outer(outer) 00176 {} 00177 00178 inline Index outer() const { return m_outer; } 00179 inline Index row() const { return m_outer; } 00180 }; 00181 00182 } // end namespace internal 00183 00184 // SparseMatrixBase functions 00185 00186 template<typename Derived> 00187 template<typename OtherDerived> 00188 const SparseDiagonalProduct<Derived,OtherDerived> 00189 SparseMatrixBase<Derived>::operator*(const DiagonalBase<OtherDerived> &other) const 00190 { 00191 return SparseDiagonalProduct<Derived,OtherDerived>(this->derived(), other.derived()); 00192 } 00193 00194 } // end namespace Eigen 00195 00196 #endif // EIGEN_SPARSE_DIAGONAL_PRODUCT_H