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 nested<MatrixType>::type MatrixTypeNested; 36 typedef typename remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned; 37 typedef MatrixType ExpressionType; 38 typedef typename MatrixType::PlainObject DenseMatrixType; 39 enum { 40 Mode = UpLo | SelfAdjoint, 41 Flags = MatrixTypeNestedCleaned::Flags & (HereditaryBits) 42 & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit)), // FIXME these flags should be preserved 43 CoeffReadCost = MatrixTypeNestedCleaned::CoeffReadCost 44 }; 45 }; 46 } 47 48 template <typename Lhs, int LhsMode, bool LhsIsVector, 49 typename Rhs, int RhsMode, bool RhsIsVector> 50 struct SelfadjointProductMatrix; 51 52 // FIXME could also be called SelfAdjointWrapper to be consistent with DiagonalWrapper ?? 53 template<typename MatrixType, unsigned int UpLo> class SelfAdjointView 54 : public TriangularBase<SelfAdjointView<MatrixType, UpLo> > 55 { 56 public: 57 58 typedef TriangularBase<SelfAdjointView> Base; 59 typedef typename internal::traits<SelfAdjointView>::MatrixTypeNested MatrixTypeNested; 60 typedef typename internal::traits<SelfAdjointView>::MatrixTypeNestedCleaned MatrixTypeNestedCleaned; 61 62 /** \brief The type of coefficients in this matrix */ 63 typedef typename internal::traits<SelfAdjointView>::Scalar Scalar; 64 65 typedef typename MatrixType::Index Index; 66 67 enum { 68 Mode = internal::traits<SelfAdjointView>::Mode 69 }; 70 typedef typename MatrixType::PlainObject PlainObject; 71 72 inline SelfAdjointView(MatrixType& matrix) : m_matrix(matrix) 73 {} 74 75 inline Index rows() const { return m_matrix.rows(); } 76 inline Index cols() const { return m_matrix.cols(); } 77 inline Index outerStride() const { return m_matrix.outerStride(); } 78 inline Index innerStride() const { return m_matrix.innerStride(); } 79 80 /** \sa MatrixBase::coeff() 81 * \warning the coordinates must fit into the referenced triangular part 82 */ 83 inline Scalar coeff(Index row, Index col) const 84 { 85 Base::check_coordinates_internal(row, col); 86 return m_matrix.coeff(row, col); 87 } 88 89 /** \sa MatrixBase::coeffRef() 90 * \warning the coordinates must fit into the referenced triangular part 91 */ 92 inline Scalar& coeffRef(Index row, Index col) 93 { 94 Base::check_coordinates_internal(row, col); 95 return m_matrix.const_cast_derived().coeffRef(row, col); 96 } 97 98 /** \internal */ 99 const MatrixTypeNestedCleaned& _expression() const { return m_matrix; } 100 101 const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; } 102 MatrixTypeNestedCleaned& nestedExpression() { return *const_cast<MatrixTypeNestedCleaned*>(&m_matrix); } 103 104 /** Efficient self-adjoint matrix times vector/matrix product */ 105 template<typename OtherDerived> 106 SelfadjointProductMatrix<MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime> 107 operator*(const MatrixBase<OtherDerived>& rhs) const 108 { 109 return SelfadjointProductMatrix 110 <MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime> 111 (m_matrix, rhs.derived()); 112 } 113 114 /** Efficient vector/matrix times self-adjoint matrix product */ 115 template<typename OtherDerived> friend 116 SelfadjointProductMatrix<OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false> 117 operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView& rhs) 118 { 119 return SelfadjointProductMatrix 120 <OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false> 121 (lhs.derived(),rhs.m_matrix); 122 } 123 124 /** Perform a symmetric rank 2 update of the selfadjoint matrix \c *this: 125 * \f$ this = this + \alpha u v^* + conj(\alpha) v u^* \f$ 126 * \returns a reference to \c *this 127 * 128 * The vectors \a u and \c v \b must be column vectors, however they can be 129 * a adjoint expression without any overhead. Only the meaningful triangular 130 * part of the matrix is updated, the rest is left unchanged. 131 * 132 * \sa rankUpdate(const MatrixBase<DerivedU>&, Scalar) 133 */ 134 template<typename DerivedU, typename DerivedV> 135 SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, Scalar alpha = Scalar(1)); 136 137 /** Perform a symmetric rank K update of the selfadjoint matrix \c *this: 138 * \f$ this = this + \alpha ( u u^* ) \f$ where \a u is a vector or matrix. 139 * 140 * \returns a reference to \c *this 141 * 142 * Note that to perform \f$ this = this + \alpha ( u^* u ) \f$ you can simply 143 * call this function with u.adjoint(). 144 * 145 * \sa rankUpdate(const MatrixBase<DerivedU>&, const MatrixBase<DerivedV>&, Scalar) 146 */ 147 template<typename DerivedU> 148 SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, Scalar alpha = Scalar(1)); 149 150 /////////// Cholesky module /////////// 151 152 const LLT<PlainObject, UpLo> llt() const; 153 const LDLT<PlainObject, UpLo> ldlt() const; 154 155 /////////// Eigenvalue module /////////// 156 157 /** Real part of #Scalar */ 158 typedef typename NumTraits<Scalar>::Real RealScalar; 159 /** Return type of eigenvalues() */ 160 typedef Matrix<RealScalar, internal::traits<MatrixType>::ColsAtCompileTime, 1> EigenvaluesReturnType; 161 162 EigenvaluesReturnType eigenvalues() const; 163 RealScalar operatorNorm() const; 164 165 #ifdef EIGEN2_SUPPORT 166 template<typename OtherDerived> 167 SelfAdjointView& operator=(const MatrixBase<OtherDerived>& other) 168 { 169 enum { 170 OtherPart = UpLo == Upper ? StrictlyLower : StrictlyUpper 171 }; 172 m_matrix.const_cast_derived().template triangularView<UpLo>() = other; 173 m_matrix.const_cast_derived().template triangularView<OtherPart>() = other.adjoint(); 174 return *this; 175 } 176 template<typename OtherMatrixType, unsigned int OtherMode> 177 SelfAdjointView& operator=(const TriangularView<OtherMatrixType, OtherMode>& other) 178 { 179 enum { 180 OtherPart = UpLo == Upper ? StrictlyLower : StrictlyUpper 181 }; 182 m_matrix.const_cast_derived().template triangularView<UpLo>() = other.toDenseMatrix(); 183 m_matrix.const_cast_derived().template triangularView<OtherPart>() = other.toDenseMatrix().adjoint(); 184 return *this; 185 } 186 #endif 187 188 protected: 189 MatrixTypeNested m_matrix; 190 }; 191 192 193 // template<typename OtherDerived, typename MatrixType, unsigned int UpLo> 194 // internal::selfadjoint_matrix_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> > 195 // operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView<MatrixType,UpLo>& rhs) 196 // { 197 // return internal::matrix_selfadjoint_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> >(lhs.derived(),rhs); 198 // } 199 200 // selfadjoint to dense matrix 201 202 namespace internal { 203 204 template<typename Derived1, typename Derived2, int UnrollCount, bool ClearOpposite> 205 struct triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Upper), UnrollCount, ClearOpposite> 206 { 207 enum { 208 col = (UnrollCount-1) / Derived1::RowsAtCompileTime, 209 row = (UnrollCount-1) % Derived1::RowsAtCompileTime 210 }; 211 212 static inline void run(Derived1 &dst, const Derived2 &src) 213 { 214 triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Upper), UnrollCount-1, ClearOpposite>::run(dst, src); 215 216 if(row == col) 217 dst.coeffRef(row, col) = real(src.coeff(row, col)); 218 else if(row < col) 219 dst.coeffRef(col, row) = conj(dst.coeffRef(row, col) = src.coeff(row, col)); 220 } 221 }; 222 223 template<typename Derived1, typename Derived2, bool ClearOpposite> 224 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Upper, 0, ClearOpposite> 225 { 226 static inline void run(Derived1 &, const Derived2 &) {} 227 }; 228 229 template<typename Derived1, typename Derived2, int UnrollCount, bool ClearOpposite> 230 struct triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Lower), UnrollCount, ClearOpposite> 231 { 232 enum { 233 col = (UnrollCount-1) / Derived1::RowsAtCompileTime, 234 row = (UnrollCount-1) % Derived1::RowsAtCompileTime 235 }; 236 237 static inline void run(Derived1 &dst, const Derived2 &src) 238 { 239 triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Lower), UnrollCount-1, ClearOpposite>::run(dst, src); 240 241 if(row == col) 242 dst.coeffRef(row, col) = real(src.coeff(row, col)); 243 else if(row > col) 244 dst.coeffRef(col, row) = conj(dst.coeffRef(row, col) = src.coeff(row, col)); 245 } 246 }; 247 248 template<typename Derived1, typename Derived2, bool ClearOpposite> 249 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Lower, 0, ClearOpposite> 250 { 251 static inline void run(Derived1 &, const Derived2 &) {} 252 }; 253 254 template<typename Derived1, typename Derived2, bool ClearOpposite> 255 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Upper, Dynamic, ClearOpposite> 256 { 257 typedef typename Derived1::Index Index; 258 static inline void run(Derived1 &dst, const Derived2 &src) 259 { 260 for(Index j = 0; j < dst.cols(); ++j) 261 { 262 for(Index i = 0; i < j; ++i) 263 { 264 dst.copyCoeff(i, j, src); 265 dst.coeffRef(j,i) = conj(dst.coeff(i,j)); 266 } 267 dst.copyCoeff(j, j, src); 268 } 269 } 270 }; 271 272 template<typename Derived1, typename Derived2, bool ClearOpposite> 273 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Lower, Dynamic, ClearOpposite> 274 { 275 static inline void run(Derived1 &dst, const Derived2 &src) 276 { 277 typedef typename Derived1::Index Index; 278 for(Index i = 0; i < dst.rows(); ++i) 279 { 280 for(Index j = 0; j < i; ++j) 281 { 282 dst.copyCoeff(i, j, src); 283 dst.coeffRef(j,i) = conj(dst.coeff(i,j)); 284 } 285 dst.copyCoeff(i, i, src); 286 } 287 } 288 }; 289 290 } // end namespace internal 291 292 /*************************************************************************** 293 * Implementation of MatrixBase methods 294 ***************************************************************************/ 295 296 template<typename Derived> 297 template<unsigned int UpLo> 298 typename MatrixBase<Derived>::template ConstSelfAdjointViewReturnType<UpLo>::Type 299 MatrixBase<Derived>::selfadjointView() const 300 { 301 return derived(); 302 } 303 304 template<typename Derived> 305 template<unsigned int UpLo> 306 typename MatrixBase<Derived>::template SelfAdjointViewReturnType<UpLo>::Type 307 MatrixBase<Derived>::selfadjointView() 308 { 309 return derived(); 310 } 311 312 } // end namespace Eigen 313 314 #endif // EIGEN_SELFADJOINTMATRIX_H 315