$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) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 // Copyright (C) 2007-2009 Benoit Jacob <jacob.benoit.1@gmail.com> 00006 // 00007 // This Source Code Form is subject to the terms of the Mozilla 00008 // Public License v. 2.0. If a copy of the MPL was not distributed 00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00010 00011 #ifndef EIGEN_DIAGONALPRODUCT_H 00012 #define EIGEN_DIAGONALPRODUCT_H 00013 00014 namespace Eigen { 00015 00016 namespace internal { 00017 template<typename MatrixType, typename DiagonalType, int ProductOrder> 00018 struct traits<DiagonalProduct<MatrixType, DiagonalType, ProductOrder> > 00019 : traits<MatrixType> 00020 { 00021 typedef typename scalar_product_traits<typename MatrixType::Scalar, typename DiagonalType::Scalar>::ReturnType Scalar; 00022 enum { 00023 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00024 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00025 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00026 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 00027 00028 _StorageOrder = MatrixType::Flags & RowMajorBit ? RowMajor : ColMajor, 00029 _ScalarAccessOnDiag = !((int(_StorageOrder) == ColMajor && int(ProductOrder) == OnTheLeft) 00030 ||(int(_StorageOrder) == RowMajor && int(ProductOrder) == OnTheRight)), 00031 _SameTypes = is_same<typename MatrixType::Scalar, typename DiagonalType::Scalar>::value, 00032 // FIXME currently we need same types, but in the future the next rule should be the one 00033 //_Vectorizable = bool(int(MatrixType::Flags)&PacketAccessBit) && ((!_PacketOnDiag) || (_SameTypes && bool(int(DiagonalType::DiagonalVectorType::Flags)&PacketAccessBit))), 00034 _Vectorizable = bool(int(MatrixType::Flags)&PacketAccessBit) && _SameTypes && (_ScalarAccessOnDiag || (bool(int(DiagonalType::DiagonalVectorType::Flags)&PacketAccessBit))), 00035 _LinearAccessMask = (RowsAtCompileTime==1 || ColsAtCompileTime==1) ? LinearAccessBit : 0, 00036 00037 Flags = ((HereditaryBits|_LinearAccessMask|AlignedBit) & (unsigned int)(MatrixType::Flags)) | (_Vectorizable ? PacketAccessBit : 0),//(int(MatrixType::Flags)&int(DiagonalType::DiagonalVectorType::Flags)&AlignedBit), 00038 CoeffReadCost = NumTraits<Scalar>::MulCost + MatrixType::CoeffReadCost + DiagonalType::DiagonalVectorType::CoeffReadCost 00039 }; 00040 }; 00041 } 00042 00043 template<typename MatrixType, typename DiagonalType, int ProductOrder> 00044 class DiagonalProduct : internal::no_assignment_operator, 00045 public MatrixBase<DiagonalProduct<MatrixType, DiagonalType, ProductOrder> > 00046 { 00047 public: 00048 00049 typedef MatrixBase<DiagonalProduct> Base; 00050 EIGEN_DENSE_PUBLIC_INTERFACE(DiagonalProduct) 00051 00052 inline DiagonalProduct(const MatrixType& matrix, const DiagonalType& diagonal) 00053 : m_matrix(matrix), m_diagonal(diagonal) 00054 { 00055 eigen_assert(diagonal.diagonal().size() == (ProductOrder == OnTheLeft ? matrix.rows() : matrix.cols())); 00056 } 00057 00058 EIGEN_STRONG_INLINE Index rows() const { return m_matrix.rows(); } 00059 EIGEN_STRONG_INLINE Index cols() const { return m_matrix.cols(); } 00060 00061 EIGEN_STRONG_INLINE const Scalar coeff(Index row, Index col) const 00062 { 00063 return m_diagonal.diagonal().coeff(ProductOrder == OnTheLeft ? row : col) * m_matrix.coeff(row, col); 00064 } 00065 00066 EIGEN_STRONG_INLINE const Scalar coeff(Index idx) const 00067 { 00068 enum { 00069 StorageOrder = int(MatrixType::Flags) & RowMajorBit ? RowMajor : ColMajor 00070 }; 00071 return coeff(int(StorageOrder)==ColMajor?idx:0,int(StorageOrder)==ColMajor?0:idx); 00072 } 00073 00074 template<int LoadMode> 00075 EIGEN_STRONG_INLINE PacketScalar packet(Index row, Index col) const 00076 { 00077 enum { 00078 StorageOrder = Flags & RowMajorBit ? RowMajor : ColMajor 00079 }; 00080 const Index indexInDiagonalVector = ProductOrder == OnTheLeft ? row : col; 00081 return packet_impl<LoadMode>(row,col,indexInDiagonalVector,typename internal::conditional< 00082 ((int(StorageOrder) == RowMajor && int(ProductOrder) == OnTheLeft) 00083 ||(int(StorageOrder) == ColMajor && int(ProductOrder) == OnTheRight)), internal::true_type, internal::false_type>::type()); 00084 } 00085 00086 template<int LoadMode> 00087 EIGEN_STRONG_INLINE PacketScalar packet(Index idx) const 00088 { 00089 enum { 00090 StorageOrder = int(MatrixType::Flags) & RowMajorBit ? RowMajor : ColMajor 00091 }; 00092 return packet<LoadMode>(int(StorageOrder)==ColMajor?idx:0,int(StorageOrder)==ColMajor?0:idx); 00093 } 00094 00095 protected: 00096 template<int LoadMode> 00097 EIGEN_STRONG_INLINE PacketScalar packet_impl(Index row, Index col, Index id, internal::true_type) const 00098 { 00099 return internal::pmul(m_matrix.template packet<LoadMode>(row, col), 00100 internal::pset1<PacketScalar>(m_diagonal.diagonal().coeff(id))); 00101 } 00102 00103 template<int LoadMode> 00104 EIGEN_STRONG_INLINE PacketScalar packet_impl(Index row, Index col, Index id, internal::false_type) const 00105 { 00106 enum { 00107 InnerSize = (MatrixType::Flags & RowMajorBit) ? MatrixType::ColsAtCompileTime : MatrixType::RowsAtCompileTime, 00108 DiagonalVectorPacketLoadMode = (LoadMode == Aligned && (((InnerSize%16) == 0) || (int(DiagonalType::DiagonalVectorType::Flags)&AlignedBit)==AlignedBit) ? Aligned : Unaligned) 00109 }; 00110 return internal::pmul(m_matrix.template packet<LoadMode>(row, col), 00111 m_diagonal.diagonal().template packet<DiagonalVectorPacketLoadMode>(id)); 00112 } 00113 00114 typename MatrixType::Nested m_matrix; 00115 typename DiagonalType::Nested m_diagonal; 00116 }; 00117 00120 template<typename Derived> 00121 template<typename DiagonalDerived> 00122 inline const DiagonalProduct<Derived, DiagonalDerived, OnTheRight> 00123 MatrixBase<Derived>::operator*(const DiagonalBase<DiagonalDerived> &a_diagonal) const 00124 { 00125 return DiagonalProduct<Derived, DiagonalDerived, OnTheRight>(derived(), a_diagonal.derived()); 00126 } 00127 00128 } // end namespace Eigen 00129 00130 #endif // EIGEN_DIAGONALPRODUCT_H