SparseDiagonalProduct.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) 2009 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_SPARSE_DIAGONAL_PRODUCT_H
00026 #define EIGEN_SPARSE_DIAGONAL_PRODUCT_H
00027 
00028 namespace Eigen { 
00029 
00030 // The product of a diagonal matrix with a sparse matrix can be easily
00031 // implemented using expression template.
00032 // We have two consider very different cases:
00033 // 1 - diag * row-major sparse
00034 //     => each inner vector <=> scalar * sparse vector product
00035 //     => so we can reuse CwiseUnaryOp::InnerIterator
00036 // 2 - diag * col-major sparse
00037 //     => each inner vector <=> densevector * sparse vector cwise product
00038 //     => again, we can reuse specialization of CwiseBinaryOp::InnerIterator
00039 //        for that particular case
00040 // The two other cases are symmetric.
00041 
00042 namespace internal {
00043 
00044 template<typename Lhs, typename Rhs>
00045 struct traits<SparseDiagonalProduct<Lhs, Rhs> >
00046 {
00047   typedef typename remove_all<Lhs>::type _Lhs;
00048   typedef typename remove_all<Rhs>::type _Rhs;
00049   typedef typename _Lhs::Scalar Scalar;
00050   typedef typename promote_index_type<typename traits<Lhs>::Index,
00051                                          typename traits<Rhs>::Index>::type Index;
00052   typedef Sparse StorageKind;
00053   typedef MatrixXpr XprKind;
00054   enum {
00055     RowsAtCompileTime = _Lhs::RowsAtCompileTime,
00056     ColsAtCompileTime = _Rhs::ColsAtCompileTime,
00057 
00058     MaxRowsAtCompileTime = _Lhs::MaxRowsAtCompileTime,
00059     MaxColsAtCompileTime = _Rhs::MaxColsAtCompileTime,
00060 
00061     SparseFlags = is_diagonal<_Lhs>::ret ? int(_Rhs::Flags) : int(_Lhs::Flags),
00062     Flags = (SparseFlags&RowMajorBit),
00063     CoeffReadCost = Dynamic
00064   };
00065 };
00066 
00067 enum {SDP_IsDiagonal, SDP_IsSparseRowMajor, SDP_IsSparseColMajor};
00068 template<typename Lhs, typename Rhs, typename SparseDiagonalProductType, int RhsMode, int LhsMode>
00069 class sparse_diagonal_product_inner_iterator_selector;
00070 
00071 } // end namespace internal
00072 
00073 template<typename Lhs, typename Rhs>
00074 class SparseDiagonalProduct
00075   : public SparseMatrixBase<SparseDiagonalProduct<Lhs,Rhs> >,
00076     internal::no_assignment_operator
00077 {
00078     typedef typename Lhs::Nested LhsNested;
00079     typedef typename Rhs::Nested RhsNested;
00080 
00081     typedef typename internal::remove_all<LhsNested>::type _LhsNested;
00082     typedef typename internal::remove_all<RhsNested>::type _RhsNested;
00083 
00084     enum {
00085       LhsMode = internal::is_diagonal<_LhsNested>::ret ? internal::SDP_IsDiagonal
00086               : (_LhsNested::Flags&RowMajorBit) ? internal::SDP_IsSparseRowMajor : internal::SDP_IsSparseColMajor,
00087       RhsMode = internal::is_diagonal<_RhsNested>::ret ? internal::SDP_IsDiagonal
00088               : (_RhsNested::Flags&RowMajorBit) ? internal::SDP_IsSparseRowMajor : internal::SDP_IsSparseColMajor
00089     };
00090 
00091   public:
00092 
00093     EIGEN_SPARSE_PUBLIC_INTERFACE(SparseDiagonalProduct)
00094 
00095     typedef internal::sparse_diagonal_product_inner_iterator_selector
00096                 <_LhsNested,_RhsNested,SparseDiagonalProduct,LhsMode,RhsMode> InnerIterator;
00097 
00098     EIGEN_STRONG_INLINE SparseDiagonalProduct(const Lhs& lhs, const Rhs& rhs)
00099       : m_lhs(lhs), m_rhs(rhs)
00100     {
00101       eigen_assert(lhs.cols() == rhs.rows() && "invalid sparse matrix * diagonal matrix product");
00102     }
00103 
00104     EIGEN_STRONG_INLINE Index rows() const { return m_lhs.rows(); }
00105     EIGEN_STRONG_INLINE Index cols() const { return m_rhs.cols(); }
00106 
00107     EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; }
00108     EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; }
00109 
00110   protected:
00111     LhsNested m_lhs;
00112     RhsNested m_rhs;
00113 };
00114 
00115 namespace internal {
00116 
00117 template<typename Lhs, typename Rhs, typename SparseDiagonalProductType>
00118 class sparse_diagonal_product_inner_iterator_selector
00119 <Lhs,Rhs,SparseDiagonalProductType,SDP_IsDiagonal,SDP_IsSparseRowMajor>
00120   : public CwiseUnaryOp<scalar_multiple_op<typename Lhs::Scalar>,const Rhs>::InnerIterator
00121 {
00122     typedef typename CwiseUnaryOp<scalar_multiple_op<typename Lhs::Scalar>,const Rhs>::InnerIterator Base;
00123     typedef typename Lhs::Index Index;
00124   public:
00125     inline sparse_diagonal_product_inner_iterator_selector(
00126               const SparseDiagonalProductType& expr, Index outer)
00127       : Base(expr.rhs()*(expr.lhs().diagonal().coeff(outer)), outer)
00128     {}
00129 };
00130 
00131 template<typename Lhs, typename Rhs, typename SparseDiagonalProductType>
00132 class sparse_diagonal_product_inner_iterator_selector
00133 <Lhs,Rhs,SparseDiagonalProductType,SDP_IsDiagonal,SDP_IsSparseColMajor>
00134   : public CwiseBinaryOp<
00135       scalar_product_op<typename Lhs::Scalar>,
00136       SparseInnerVectorSet<Rhs,1>,
00137       typename Lhs::DiagonalVectorType>::InnerIterator
00138 {
00139     typedef typename CwiseBinaryOp<
00140       scalar_product_op<typename Lhs::Scalar>,
00141       SparseInnerVectorSet<Rhs,1>,
00142       typename Lhs::DiagonalVectorType>::InnerIterator Base;
00143     typedef typename Lhs::Index Index;
00144   public:
00145     inline sparse_diagonal_product_inner_iterator_selector(
00146               const SparseDiagonalProductType& expr, Index outer)
00147       : Base(expr.rhs().innerVector(outer) .cwiseProduct(expr.lhs().diagonal()), 0)
00148     {}
00149 };
00150 
00151 template<typename Lhs, typename Rhs, typename SparseDiagonalProductType>
00152 class sparse_diagonal_product_inner_iterator_selector
00153 <Lhs,Rhs,SparseDiagonalProductType,SDP_IsSparseColMajor,SDP_IsDiagonal>
00154   : public CwiseUnaryOp<scalar_multiple_op<typename Rhs::Scalar>,const Lhs>::InnerIterator
00155 {
00156     typedef typename CwiseUnaryOp<scalar_multiple_op<typename Rhs::Scalar>,const Lhs>::InnerIterator Base;
00157     typedef typename Lhs::Index Index;
00158   public:
00159     inline sparse_diagonal_product_inner_iterator_selector(
00160               const SparseDiagonalProductType& expr, Index outer)
00161       : Base(expr.lhs()*expr.rhs().diagonal().coeff(outer), outer)
00162     {}
00163 };
00164 
00165 template<typename Lhs, typename Rhs, typename SparseDiagonalProductType>
00166 class sparse_diagonal_product_inner_iterator_selector
00167 <Lhs,Rhs,SparseDiagonalProductType,SDP_IsSparseRowMajor,SDP_IsDiagonal>
00168   : public CwiseBinaryOp<
00169       scalar_product_op<typename Rhs::Scalar>,
00170       SparseInnerVectorSet<Lhs,1>,
00171       Transpose<const typename Rhs::DiagonalVectorType> >::InnerIterator
00172 {
00173     typedef typename CwiseBinaryOp<
00174       scalar_product_op<typename Rhs::Scalar>,
00175       SparseInnerVectorSet<Lhs,1>,
00176       Transpose<const typename Rhs::DiagonalVectorType> >::InnerIterator Base;
00177     typedef typename Lhs::Index Index;
00178   public:
00179     inline sparse_diagonal_product_inner_iterator_selector(
00180               const SparseDiagonalProductType& expr, Index outer)
00181       : Base(expr.lhs().innerVector(outer) .cwiseProduct(expr.rhs().diagonal().transpose()), 0)
00182     {}
00183 };
00184 
00185 } // end namespace internal
00186 
00187 // SparseMatrixBase functions
00188 
00189 template<typename Derived>
00190 template<typename OtherDerived>
00191 const SparseDiagonalProduct<Derived,OtherDerived>
00192 SparseMatrixBase<Derived>::operator*(const DiagonalBase<OtherDerived> &other) const
00193 {
00194   return SparseDiagonalProduct<Derived,OtherDerived>(this->derived(), other.derived());
00195 }
00196 
00197 } // end namespace Eigen
00198 
00199 #endif // EIGEN_SPARSE_DIAGONAL_PRODUCT_H