![]() |
Eigen
3.3.3
|
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_SELFADJOINTMATRIX_H 00011 #define EIGEN_SELFADJOINTMATRIX_H 00012 00013 namespace Eigen { 00014 00031 namespace internal { 00032 template<typename MatrixType, unsigned int UpLo> 00033 struct traits<SelfAdjointView<MatrixType, UpLo> > : traits<MatrixType> 00034 { 00035 typedef typename ref_selector<MatrixType>::non_const_type MatrixTypeNested; 00036 typedef typename remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned; 00037 typedef MatrixType ExpressionType; 00038 typedef typename MatrixType::PlainObject FullMatrixType; 00039 enum { 00040 Mode = UpLo | SelfAdjoint, 00041 FlagsLvalueBit = is_lvalue<MatrixType>::value ? LvalueBit : 0, 00042 Flags = MatrixTypeNestedCleaned::Flags & (HereditaryBits|FlagsLvalueBit) 00043 & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit)) // FIXME these flags should be preserved 00044 }; 00045 }; 00046 } 00047 00048 00049 template<typename _MatrixType, unsigned int UpLo> class SelfAdjointView 00050 : public TriangularBase<SelfAdjointView<_MatrixType, UpLo> > 00051 { 00052 public: 00053 00054 typedef _MatrixType MatrixType; 00055 typedef TriangularBase<SelfAdjointView> Base; 00056 typedef typename internal::traits<SelfAdjointView>::MatrixTypeNested MatrixTypeNested; 00057 typedef typename internal::traits<SelfAdjointView>::MatrixTypeNestedCleaned MatrixTypeNestedCleaned; 00058 typedef MatrixTypeNestedCleaned NestedExpression; 00059 00061 typedef typename internal::traits<SelfAdjointView>::Scalar Scalar; 00062 typedef typename MatrixType::StorageIndex StorageIndex; 00063 typedef typename internal::remove_all<typename MatrixType::ConjugateReturnType>::type MatrixConjugateReturnType; 00064 00065 enum { 00066 Mode = internal::traits<SelfAdjointView>::Mode, 00067 Flags = internal::traits<SelfAdjointView>::Flags, 00068 TransposeMode = ((Mode & Upper) ? Lower : 0) | ((Mode & Lower) ? Upper : 0) 00069 }; 00070 typedef typename MatrixType::PlainObject PlainObject; 00071 00072 EIGEN_DEVICE_FUNC 00073 explicit inline SelfAdjointView(MatrixType& matrix) : m_matrix(matrix) 00074 {} 00075 00076 EIGEN_DEVICE_FUNC 00077 inline Index rows() const { return m_matrix.rows(); } 00078 EIGEN_DEVICE_FUNC 00079 inline Index cols() const { return m_matrix.cols(); } 00080 EIGEN_DEVICE_FUNC 00081 inline Index outerStride() const { return m_matrix.outerStride(); } 00082 EIGEN_DEVICE_FUNC 00083 inline Index innerStride() const { return m_matrix.innerStride(); } 00084 00088 EIGEN_DEVICE_FUNC 00089 inline Scalar coeff(Index row, Index col) const 00090 { 00091 Base::check_coordinates_internal(row, col); 00092 return m_matrix.coeff(row, col); 00093 } 00094 00098 EIGEN_DEVICE_FUNC 00099 inline Scalar& coeffRef(Index row, Index col) 00100 { 00101 EIGEN_STATIC_ASSERT_LVALUE(SelfAdjointView); 00102 Base::check_coordinates_internal(row, col); 00103 return m_matrix.coeffRef(row, col); 00104 } 00105 00107 EIGEN_DEVICE_FUNC 00108 const MatrixTypeNestedCleaned& _expression() const { return m_matrix; } 00109 00110 EIGEN_DEVICE_FUNC 00111 const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; } 00112 EIGEN_DEVICE_FUNC 00113 MatrixTypeNestedCleaned& nestedExpression() { return m_matrix; } 00114 00116 template<typename OtherDerived> 00117 EIGEN_DEVICE_FUNC 00118 const Product<SelfAdjointView,OtherDerived> 00119 operator*(const MatrixBase<OtherDerived>& rhs) const 00120 { 00121 return Product<SelfAdjointView,OtherDerived>(*this, rhs.derived()); 00122 } 00123 00125 template<typename OtherDerived> friend 00126 EIGEN_DEVICE_FUNC 00127 const Product<OtherDerived,SelfAdjointView> 00128 operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView& rhs) 00129 { 00130 return Product<OtherDerived,SelfAdjointView>(lhs.derived(),rhs); 00131 } 00132 00133 friend EIGEN_DEVICE_FUNC 00134 const SelfAdjointView<const EIGEN_SCALAR_BINARYOP_EXPR_RETURN_TYPE(Scalar,MatrixType,product),UpLo> 00135 operator*(const Scalar& s, const SelfAdjointView& mat) 00136 { 00137 return (s*mat.nestedExpression()).template selfadjointView<UpLo>(); 00138 } 00139 00150 template<typename DerivedU, typename DerivedV> 00151 EIGEN_DEVICE_FUNC 00152 SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, const Scalar& alpha = Scalar(1)); 00153 00164 template<typename DerivedU> 00165 EIGEN_DEVICE_FUNC 00166 SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const Scalar& alpha = Scalar(1)); 00167 00178 template<unsigned int TriMode> 00179 EIGEN_DEVICE_FUNC 00180 typename internal::conditional<(TriMode&(Upper|Lower))==(UpLo&(Upper|Lower)), 00181 TriangularView<MatrixType,TriMode>, 00182 TriangularView<typename MatrixType::AdjointReturnType,TriMode> >::type 00183 triangularView() const 00184 { 00185 typename internal::conditional<(TriMode&(Upper|Lower))==(UpLo&(Upper|Lower)), MatrixType&, typename MatrixType::ConstTransposeReturnType>::type tmp1(m_matrix); 00186 typename internal::conditional<(TriMode&(Upper|Lower))==(UpLo&(Upper|Lower)), MatrixType&, typename MatrixType::AdjointReturnType>::type tmp2(tmp1); 00187 return typename internal::conditional<(TriMode&(Upper|Lower))==(UpLo&(Upper|Lower)), 00188 TriangularView<MatrixType,TriMode>, 00189 TriangularView<typename MatrixType::AdjointReturnType,TriMode> >::type(tmp2); 00190 } 00191 00192 typedef SelfAdjointView<const MatrixConjugateReturnType,Mode> ConjugateReturnType; 00194 EIGEN_DEVICE_FUNC 00195 inline const ConjugateReturnType conjugate() const 00196 { return ConjugateReturnType(m_matrix.conjugate()); } 00197 00198 typedef SelfAdjointView<const typename MatrixType::AdjointReturnType,TransposeMode> AdjointReturnType; 00200 EIGEN_DEVICE_FUNC 00201 inline const AdjointReturnType adjoint() const 00202 { return AdjointReturnType(m_matrix.adjoint()); } 00203 00204 typedef SelfAdjointView<typename MatrixType::TransposeReturnType,TransposeMode> TransposeReturnType; 00206 EIGEN_DEVICE_FUNC 00207 inline TransposeReturnType transpose() 00208 { 00209 EIGEN_STATIC_ASSERT_LVALUE(MatrixType) 00210 typename MatrixType::TransposeReturnType tmp(m_matrix); 00211 return TransposeReturnType(tmp); 00212 } 00213 00214 typedef SelfAdjointView<const typename MatrixType::ConstTransposeReturnType,TransposeMode> ConstTransposeReturnType; 00216 EIGEN_DEVICE_FUNC 00217 inline const ConstTransposeReturnType transpose() const 00218 { 00219 return ConstTransposeReturnType(m_matrix.transpose()); 00220 } 00221 00227 EIGEN_DEVICE_FUNC 00228 typename MatrixType::ConstDiagonalReturnType diagonal() const 00229 { 00230 return typename MatrixType::ConstDiagonalReturnType(m_matrix); 00231 } 00232 00234 00235 const LLT<PlainObject, UpLo> llt() const; 00236 const LDLT<PlainObject, UpLo> ldlt() const; 00237 00239 00241 typedef typename NumTraits<Scalar>::Real RealScalar; 00243 typedef Matrix<RealScalar, internal::traits<MatrixType>::ColsAtCompileTime, 1> EigenvaluesReturnType; 00244 00245 EIGEN_DEVICE_FUNC 00246 EigenvaluesReturnType eigenvalues() const; 00247 EIGEN_DEVICE_FUNC 00248 RealScalar operatorNorm() const; 00249 00250 protected: 00251 MatrixTypeNested m_matrix; 00252 }; 00253 00254 00255 // template<typename OtherDerived, typename MatrixType, unsigned int UpLo> 00256 // internal::selfadjoint_matrix_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> > 00257 // operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView<MatrixType,UpLo>& rhs) 00258 // { 00259 // return internal::matrix_selfadjoint_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> >(lhs.derived(),rhs); 00260 // } 00261 00262 // selfadjoint to dense matrix 00263 00264 namespace internal { 00265 00266 // TODO currently a selfadjoint expression has the form SelfAdjointView<.,.> 00267 // in the future selfadjoint-ness should be defined by the expression traits 00268 // such that Transpose<SelfAdjointView<.,.> > is valid. (currently TriangularBase::transpose() is overloaded to make it work) 00269 template<typename MatrixType, unsigned int Mode> 00270 struct evaluator_traits<SelfAdjointView<MatrixType,Mode> > 00271 { 00272 typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind; 00273 typedef SelfAdjointShape Shape; 00274 }; 00275 00276 template<int UpLo, int SetOpposite, typename DstEvaluatorTypeT, typename SrcEvaluatorTypeT, typename Functor, int Version> 00277 class triangular_dense_assignment_kernel<UpLo,SelfAdjoint,SetOpposite,DstEvaluatorTypeT,SrcEvaluatorTypeT,Functor,Version> 00278 : public generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version> 00279 { 00280 protected: 00281 typedef generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version> Base; 00282 typedef typename Base::DstXprType DstXprType; 00283 typedef typename Base::SrcXprType SrcXprType; 00284 using Base::m_dst; 00285 using Base::m_src; 00286 using Base::m_functor; 00287 public: 00288 00289 typedef typename Base::DstEvaluatorType DstEvaluatorType; 00290 typedef typename Base::SrcEvaluatorType SrcEvaluatorType; 00291 typedef typename Base::Scalar Scalar; 00292 typedef typename Base::AssignmentTraits AssignmentTraits; 00293 00294 00295 EIGEN_DEVICE_FUNC triangular_dense_assignment_kernel(DstEvaluatorType &dst, const SrcEvaluatorType &src, const Functor &func, DstXprType& dstExpr) 00296 : Base(dst, src, func, dstExpr) 00297 {} 00298 00299 EIGEN_DEVICE_FUNC void assignCoeff(Index row, Index col) 00300 { 00301 eigen_internal_assert(row!=col); 00302 Scalar tmp = m_src.coeff(row,col); 00303 m_functor.assignCoeff(m_dst.coeffRef(row,col), tmp); 00304 m_functor.assignCoeff(m_dst.coeffRef(col,row), numext::conj(tmp)); 00305 } 00306 00307 EIGEN_DEVICE_FUNC void assignDiagonalCoeff(Index id) 00308 { 00309 Base::assignCoeff(id,id); 00310 } 00311 00312 EIGEN_DEVICE_FUNC void assignOppositeCoeff(Index, Index) 00313 { eigen_internal_assert(false && "should never be called"); } 00314 }; 00315 00316 } // end namespace internal 00317 00318 /*************************************************************************** 00319 * Implementation of MatrixBase methods 00320 ***************************************************************************/ 00321 00323 template<typename Derived> 00324 template<unsigned int UpLo> 00325 typename MatrixBase<Derived>::template ConstSelfAdjointViewReturnType<UpLo>::Type 00326 MatrixBase<Derived>::selfadjointView() const 00327 { 00328 return typename ConstSelfAdjointViewReturnType<UpLo>::Type(derived()); 00329 } 00330 00340 template<typename Derived> 00341 template<unsigned int UpLo> 00342 typename MatrixBase<Derived>::template SelfAdjointViewReturnType<UpLo>::Type 00343 MatrixBase<Derived>::selfadjointView() 00344 { 00345 return typename SelfAdjointViewReturnType<UpLo>::Type(derived()); 00346 } 00347 00348 } // end namespace Eigen 00349 00350 #endif // EIGEN_SELFADJOINTMATRIX_H