1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_SELFADJOINTMATRIX_H 11 #define EIGEN_SELFADJOINTMATRIX_H 12 13 namespace Eigen { 14 15 /** \class SelfAdjointView 16 * \ingroup Core_Module 17 * 18 * 19 * \brief Expression of a selfadjoint matrix from a triangular part of a dense matrix 20 * 21 * \param MatrixType the type of the dense matrix storing the coefficients 22 * \param TriangularPart can be either \c #Lower or \c #Upper 23 * 24 * This class is an expression of a sefladjoint matrix from a triangular part of a matrix 25 * with given dense storage of the coefficients. It is the return type of MatrixBase::selfadjointView() 26 * and most of the time this is the only way that it is used. 27 * 28 * \sa class TriangularBase, MatrixBase::selfadjointView() 29 */ 30 31 namespace internal { 32 template<typename MatrixType, unsigned int UpLo> 33 struct traits<SelfAdjointView<MatrixType, UpLo> > : traits<MatrixType> 34 { 35 typedef typename ref_selector<MatrixType>::non_const_type MatrixTypeNested; 36 typedef typename remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned; 37 typedef MatrixType ExpressionType; 38 typedef typename MatrixType::PlainObject FullMatrixType; 39 enum { 40 Mode = UpLo | SelfAdjoint, 41 FlagsLvalueBit = is_lvalue<MatrixType>::value ? LvalueBit : 0, 42 Flags = MatrixTypeNestedCleaned::Flags & (HereditaryBits|FlagsLvalueBit) 43 & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit)) // FIXME these flags should be preserved 44 }; 45 }; 46 } 47 48 49 template<typename _MatrixType, unsigned int UpLo> class SelfAdjointView 50 : public TriangularBase<SelfAdjointView<_MatrixType, UpLo> > 51 { 52 public: 53 54 typedef _MatrixType MatrixType; 55 typedef TriangularBase<SelfAdjointView> Base; 56 typedef typename internal::traits<SelfAdjointView>::MatrixTypeNested MatrixTypeNested; 57 typedef typename internal::traits<SelfAdjointView>::MatrixTypeNestedCleaned MatrixTypeNestedCleaned; 58 typedef MatrixTypeNestedCleaned NestedExpression; 59 60 /** \brief The type of coefficients in this matrix */ 61 typedef typename internal::traits<SelfAdjointView>::Scalar Scalar; 62 typedef typename MatrixType::StorageIndex StorageIndex; 63 typedef typename internal::remove_all<typename MatrixType::ConjugateReturnType>::type MatrixConjugateReturnType; 64 65 enum { 66 Mode = internal::traits<SelfAdjointView>::Mode, 67 Flags = internal::traits<SelfAdjointView>::Flags, 68 TransposeMode = ((Mode & Upper) ? Lower : 0) | ((Mode & Lower) ? Upper : 0) 69 }; 70 typedef typename MatrixType::PlainObject PlainObject; 71 72 EIGEN_DEVICE_FUNC 73 explicit inline SelfAdjointView(MatrixType& matrix) : m_matrix(matrix) 74 { 75 EIGEN_STATIC_ASSERT(UpLo==Lower || UpLo==Upper,SELFADJOINTVIEW_ACCEPTS_UPPER_AND_LOWER_MODE_ONLY); 76 } 77 78 EIGEN_DEVICE_FUNC 79 inline Index rows() const { return m_matrix.rows(); } 80 EIGEN_DEVICE_FUNC 81 inline Index cols() const { return m_matrix.cols(); } 82 EIGEN_DEVICE_FUNC 83 inline Index outerStride() const { return m_matrix.outerStride(); } 84 EIGEN_DEVICE_FUNC 85 inline Index innerStride() const { return m_matrix.innerStride(); } 86 87 /** \sa MatrixBase::coeff() 88 * \warning the coordinates must fit into the referenced triangular part 89 */ 90 EIGEN_DEVICE_FUNC 91 inline Scalar coeff(Index row, Index col) const 92 { 93 Base::check_coordinates_internal(row, col); 94 return m_matrix.coeff(row, col); 95 } 96 97 /** \sa MatrixBase::coeffRef() 98 * \warning the coordinates must fit into the referenced triangular part 99 */ 100 EIGEN_DEVICE_FUNC 101 inline Scalar& coeffRef(Index row, Index col) 102 { 103 EIGEN_STATIC_ASSERT_LVALUE(SelfAdjointView); 104 Base::check_coordinates_internal(row, col); 105 return m_matrix.coeffRef(row, col); 106 } 107 108 /** \internal */ 109 EIGEN_DEVICE_FUNC 110 const MatrixTypeNestedCleaned& _expression() const { return m_matrix; } 111 112 EIGEN_DEVICE_FUNC 113 const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; } 114 EIGEN_DEVICE_FUNC 115 MatrixTypeNestedCleaned& nestedExpression() { return m_matrix; } 116 117 /** Efficient triangular matrix times vector/matrix product */ 118 template<typename OtherDerived> 119 EIGEN_DEVICE_FUNC 120 const Product<SelfAdjointView,OtherDerived> 121 operator*(const MatrixBase<OtherDerived>& rhs) const 122 { 123 return Product<SelfAdjointView,OtherDerived>(*this, rhs.derived()); 124 } 125 126 /** Efficient vector/matrix times triangular matrix product */ 127 template<typename OtherDerived> friend 128 EIGEN_DEVICE_FUNC 129 const Product<OtherDerived,SelfAdjointView> 130 operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView& rhs) 131 { 132 return Product<OtherDerived,SelfAdjointView>(lhs.derived(),rhs); 133 } 134 135 friend EIGEN_DEVICE_FUNC 136 const SelfAdjointView<const EIGEN_SCALAR_BINARYOP_EXPR_RETURN_TYPE(Scalar,MatrixType,product),UpLo> 137 operator*(const Scalar& s, const SelfAdjointView& mat) 138 { 139 return (s*mat.nestedExpression()).template selfadjointView<UpLo>(); 140 } 141 142 /** Perform a symmetric rank 2 update of the selfadjoint matrix \c *this: 143 * \f$ this = this + \alpha u v^* + conj(\alpha) v u^* \f$ 144 * \returns a reference to \c *this 145 * 146 * The vectors \a u and \c v \b must be column vectors, however they can be 147 * a adjoint expression without any overhead. Only the meaningful triangular 148 * part of the matrix is updated, the rest is left unchanged. 149 * 150 * \sa rankUpdate(const MatrixBase<DerivedU>&, Scalar) 151 */ 152 template<typename DerivedU, typename DerivedV> 153 EIGEN_DEVICE_FUNC 154 SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, const Scalar& alpha = Scalar(1)); 155 156 /** Perform a symmetric rank K update of the selfadjoint matrix \c *this: 157 * \f$ this = this + \alpha ( u u^* ) \f$ where \a u is a vector or matrix. 158 * 159 * \returns a reference to \c *this 160 * 161 * Note that to perform \f$ this = this + \alpha ( u^* u ) \f$ you can simply 162 * call this function with u.adjoint(). 163 * 164 * \sa rankUpdate(const MatrixBase<DerivedU>&, const MatrixBase<DerivedV>&, Scalar) 165 */ 166 template<typename DerivedU> 167 EIGEN_DEVICE_FUNC 168 SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const Scalar& alpha = Scalar(1)); 169 170 /** \returns an expression of a triangular view extracted from the current selfadjoint view of a given triangular part 171 * 172 * The parameter \a TriMode can have the following values: \c #Upper, \c #StrictlyUpper, \c #UnitUpper, 173 * \c #Lower, \c #StrictlyLower, \c #UnitLower. 174 * 175 * If \c TriMode references the same triangular part than \c *this, then this method simply return a \c TriangularView of the nested expression, 176 * otherwise, the nested expression is first transposed, thus returning a \c TriangularView<Transpose<MatrixType>> object. 177 * 178 * \sa MatrixBase::triangularView(), class TriangularView 179 */ 180 template<unsigned int TriMode> 181 EIGEN_DEVICE_FUNC 182 typename internal::conditional<(TriMode&(Upper|Lower))==(UpLo&(Upper|Lower)), 183 TriangularView<MatrixType,TriMode>, 184 TriangularView<typename MatrixType::AdjointReturnType,TriMode> >::type 185 triangularView() const 186 { 187 typename internal::conditional<(TriMode&(Upper|Lower))==(UpLo&(Upper|Lower)), MatrixType&, typename MatrixType::ConstTransposeReturnType>::type tmp1(m_matrix); 188 typename internal::conditional<(TriMode&(Upper|Lower))==(UpLo&(Upper|Lower)), MatrixType&, typename MatrixType::AdjointReturnType>::type tmp2(tmp1); 189 return typename internal::conditional<(TriMode&(Upper|Lower))==(UpLo&(Upper|Lower)), 190 TriangularView<MatrixType,TriMode>, 191 TriangularView<typename MatrixType::AdjointReturnType,TriMode> >::type(tmp2); 192 } 193 194 typedef SelfAdjointView<const MatrixConjugateReturnType,UpLo> ConjugateReturnType; 195 /** \sa MatrixBase::conjugate() const */ 196 EIGEN_DEVICE_FUNC 197 inline const ConjugateReturnType conjugate() const 198 { return ConjugateReturnType(m_matrix.conjugate()); } 199 200 typedef SelfAdjointView<const typename MatrixType::AdjointReturnType,TransposeMode> AdjointReturnType; 201 /** \sa MatrixBase::adjoint() const */ 202 EIGEN_DEVICE_FUNC 203 inline const AdjointReturnType adjoint() const 204 { return AdjointReturnType(m_matrix.adjoint()); } 205 206 typedef SelfAdjointView<typename MatrixType::TransposeReturnType,TransposeMode> TransposeReturnType; 207 /** \sa MatrixBase::transpose() */ 208 EIGEN_DEVICE_FUNC 209 inline TransposeReturnType transpose() 210 { 211 EIGEN_STATIC_ASSERT_LVALUE(MatrixType) 212 typename MatrixType::TransposeReturnType tmp(m_matrix); 213 return TransposeReturnType(tmp); 214 } 215 216 typedef SelfAdjointView<const typename MatrixType::ConstTransposeReturnType,TransposeMode> ConstTransposeReturnType; 217 /** \sa MatrixBase::transpose() const */ 218 EIGEN_DEVICE_FUNC 219 inline const ConstTransposeReturnType transpose() const 220 { 221 return ConstTransposeReturnType(m_matrix.transpose()); 222 } 223 224 /** \returns a const expression of the main diagonal of the matrix \c *this 225 * 226 * This method simply returns the diagonal of the nested expression, thus by-passing the SelfAdjointView decorator. 227 * 228 * \sa MatrixBase::diagonal(), class Diagonal */ 229 EIGEN_DEVICE_FUNC 230 typename MatrixType::ConstDiagonalReturnType diagonal() const 231 { 232 return typename MatrixType::ConstDiagonalReturnType(m_matrix); 233 } 234 235 /////////// Cholesky module /////////// 236 237 const LLT<PlainObject, UpLo> llt() const; 238 const LDLT<PlainObject, UpLo> ldlt() const; 239 240 /////////// Eigenvalue module /////////// 241 242 /** Real part of #Scalar */ 243 typedef typename NumTraits<Scalar>::Real RealScalar; 244 /** Return type of eigenvalues() */ 245 typedef Matrix<RealScalar, internal::traits<MatrixType>::ColsAtCompileTime, 1> EigenvaluesReturnType; 246 247 EIGEN_DEVICE_FUNC 248 EigenvaluesReturnType eigenvalues() const; 249 EIGEN_DEVICE_FUNC 250 RealScalar operatorNorm() const; 251 252 protected: 253 MatrixTypeNested m_matrix; 254 }; 255 256 257 // template<typename OtherDerived, typename MatrixType, unsigned int UpLo> 258 // internal::selfadjoint_matrix_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> > 259 // operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView<MatrixType,UpLo>& rhs) 260 // { 261 // return internal::matrix_selfadjoint_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> >(lhs.derived(),rhs); 262 // } 263 264 // selfadjoint to dense matrix 265 266 namespace internal { 267 268 // TODO currently a selfadjoint expression has the form SelfAdjointView<.,.> 269 // in the future selfadjoint-ness should be defined by the expression traits 270 // such that Transpose<SelfAdjointView<.,.> > is valid. (currently TriangularBase::transpose() is overloaded to make it work) 271 template<typename MatrixType, unsigned int Mode> 272 struct evaluator_traits<SelfAdjointView<MatrixType,Mode> > 273 { 274 typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind; 275 typedef SelfAdjointShape Shape; 276 }; 277 278 template<int UpLo, int SetOpposite, typename DstEvaluatorTypeT, typename SrcEvaluatorTypeT, typename Functor, int Version> 279 class triangular_dense_assignment_kernel<UpLo,SelfAdjoint,SetOpposite,DstEvaluatorTypeT,SrcEvaluatorTypeT,Functor,Version> 280 : public generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version> 281 { 282 protected: 283 typedef generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version> Base; 284 typedef typename Base::DstXprType DstXprType; 285 typedef typename Base::SrcXprType SrcXprType; 286 using Base::m_dst; 287 using Base::m_src; 288 using Base::m_functor; 289 public: 290 291 typedef typename Base::DstEvaluatorType DstEvaluatorType; 292 typedef typename Base::SrcEvaluatorType SrcEvaluatorType; 293 typedef typename Base::Scalar Scalar; 294 typedef typename Base::AssignmentTraits AssignmentTraits; 295 296 297 EIGEN_DEVICE_FUNC triangular_dense_assignment_kernel(DstEvaluatorType &dst, const SrcEvaluatorType &src, const Functor &func, DstXprType& dstExpr) 298 : Base(dst, src, func, dstExpr) 299 {} 300 301 EIGEN_DEVICE_FUNC void assignCoeff(Index row, Index col) 302 { 303 eigen_internal_assert(row!=col); 304 Scalar tmp = m_src.coeff(row,col); 305 m_functor.assignCoeff(m_dst.coeffRef(row,col), tmp); 306 m_functor.assignCoeff(m_dst.coeffRef(col,row), numext::conj(tmp)); 307 } 308 309 EIGEN_DEVICE_FUNC void assignDiagonalCoeff(Index id) 310 { 311 Base::assignCoeff(id,id); 312 } 313 314 EIGEN_DEVICE_FUNC void assignOppositeCoeff(Index, Index) 315 { eigen_internal_assert(false && "should never be called"); } 316 }; 317 318 } // end namespace internal 319 320 /*************************************************************************** 321 * Implementation of MatrixBase methods 322 ***************************************************************************/ 323 324 /** This is the const version of MatrixBase::selfadjointView() */ 325 template<typename Derived> 326 template<unsigned int UpLo> 327 typename MatrixBase<Derived>::template ConstSelfAdjointViewReturnType<UpLo>::Type 328 MatrixBase<Derived>::selfadjointView() const 329 { 330 return typename ConstSelfAdjointViewReturnType<UpLo>::Type(derived()); 331 } 332 333 /** \returns an expression of a symmetric/self-adjoint view extracted from the upper or lower triangular part of the current matrix 334 * 335 * The parameter \a UpLo can be either \c #Upper or \c #Lower 336 * 337 * Example: \include MatrixBase_selfadjointView.cpp 338 * Output: \verbinclude MatrixBase_selfadjointView.out 339 * 340 * \sa class SelfAdjointView 341 */ 342 template<typename Derived> 343 template<unsigned int UpLo> 344 typename MatrixBase<Derived>::template SelfAdjointViewReturnType<UpLo>::Type 345 MatrixBase<Derived>::selfadjointView() 346 { 347 return typename SelfAdjointViewReturnType<UpLo>::Type(derived()); 348 } 349 350 } // end namespace Eigen 351 352 #endif // EIGEN_SELFADJOINTMATRIX_H 353