$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) 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 nested<MatrixType>::type MatrixTypeNested; 00036 typedef typename remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned; 00037 typedef MatrixType ExpressionType; 00038 typedef typename MatrixType::PlainObject DenseMatrixType; 00039 enum { 00040 Mode = UpLo | SelfAdjoint, 00041 Flags = MatrixTypeNestedCleaned::Flags & (HereditaryBits) 00042 & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit)), // FIXME these flags should be preserved 00043 CoeffReadCost = MatrixTypeNestedCleaned::CoeffReadCost 00044 }; 00045 }; 00046 } 00047 00048 template <typename Lhs, int LhsMode, bool LhsIsVector, 00049 typename Rhs, int RhsMode, bool RhsIsVector> 00050 struct SelfadjointProductMatrix; 00051 00052 // FIXME could also be called SelfAdjointWrapper to be consistent with DiagonalWrapper ?? 00053 template<typename MatrixType, unsigned int UpLo> class SelfAdjointView 00054 : public TriangularBase<SelfAdjointView<MatrixType, UpLo> > 00055 { 00056 public: 00057 00058 typedef TriangularBase<SelfAdjointView> Base; 00059 typedef typename internal::traits<SelfAdjointView>::MatrixTypeNested MatrixTypeNested; 00060 typedef typename internal::traits<SelfAdjointView>::MatrixTypeNestedCleaned MatrixTypeNestedCleaned; 00061 00063 typedef typename internal::traits<SelfAdjointView>::Scalar Scalar; 00064 00065 typedef typename MatrixType::Index Index; 00066 00067 enum { 00068 Mode = internal::traits<SelfAdjointView>::Mode 00069 }; 00070 typedef typename MatrixType::PlainObject PlainObject; 00071 00072 inline SelfAdjointView(MatrixType& matrix) : m_matrix(matrix) 00073 {} 00074 00075 inline Index rows() const { return m_matrix.rows(); } 00076 inline Index cols() const { return m_matrix.cols(); } 00077 inline Index outerStride() const { return m_matrix.outerStride(); } 00078 inline Index innerStride() const { return m_matrix.innerStride(); } 00079 00083 inline Scalar coeff(Index row, Index col) const 00084 { 00085 Base::check_coordinates_internal(row, col); 00086 return m_matrix.coeff(row, col); 00087 } 00088 00092 inline Scalar& coeffRef(Index row, Index col) 00093 { 00094 Base::check_coordinates_internal(row, col); 00095 return m_matrix.const_cast_derived().coeffRef(row, col); 00096 } 00097 00099 const MatrixTypeNestedCleaned& _expression() const { return m_matrix; } 00100 00101 const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; } 00102 MatrixTypeNestedCleaned& nestedExpression() { return *const_cast<MatrixTypeNestedCleaned*>(&m_matrix); } 00103 00105 template<typename OtherDerived> 00106 SelfadjointProductMatrix<MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime> 00107 operator*(const MatrixBase<OtherDerived>& rhs) const 00108 { 00109 return SelfadjointProductMatrix 00110 <MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime> 00111 (m_matrix, rhs.derived()); 00112 } 00113 00115 template<typename OtherDerived> friend 00116 SelfadjointProductMatrix<OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false> 00117 operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView& rhs) 00118 { 00119 return SelfadjointProductMatrix 00120 <OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false> 00121 (lhs.derived(),rhs.m_matrix); 00122 } 00123 00134 template<typename DerivedU, typename DerivedV> 00135 SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, const Scalar& alpha = Scalar(1)); 00136 00147 template<typename DerivedU> 00148 SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const Scalar& alpha = Scalar(1)); 00149 00151 00152 const LLT<PlainObject, UpLo> llt() const; 00153 const LDLT<PlainObject, UpLo> ldlt() const; 00154 00156 00158 typedef typename NumTraits<Scalar>::Real RealScalar; 00160 typedef Matrix<RealScalar, internal::traits<MatrixType>::ColsAtCompileTime, 1> EigenvaluesReturnType; 00161 00162 EigenvaluesReturnType eigenvalues() const; 00163 RealScalar operatorNorm() const; 00164 00165 #ifdef EIGEN2_SUPPORT 00166 template<typename OtherDerived> 00167 SelfAdjointView& operator=(const MatrixBase<OtherDerived>& other) 00168 { 00169 enum { 00170 OtherPart = UpLo == Upper ? StrictlyLower : StrictlyUpper 00171 }; 00172 m_matrix.const_cast_derived().template triangularView<UpLo>() = other; 00173 m_matrix.const_cast_derived().template triangularView<OtherPart>() = other.adjoint(); 00174 return *this; 00175 } 00176 template<typename OtherMatrixType, unsigned int OtherMode> 00177 SelfAdjointView& operator=(const TriangularView<OtherMatrixType, OtherMode>& other) 00178 { 00179 enum { 00180 OtherPart = UpLo == Upper ? StrictlyLower : StrictlyUpper 00181 }; 00182 m_matrix.const_cast_derived().template triangularView<UpLo>() = other.toDenseMatrix(); 00183 m_matrix.const_cast_derived().template triangularView<OtherPart>() = other.toDenseMatrix().adjoint(); 00184 return *this; 00185 } 00186 #endif 00187 00188 protected: 00189 MatrixTypeNested m_matrix; 00190 }; 00191 00192 00193 // template<typename OtherDerived, typename MatrixType, unsigned int UpLo> 00194 // internal::selfadjoint_matrix_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> > 00195 // operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView<MatrixType,UpLo>& rhs) 00196 // { 00197 // return internal::matrix_selfadjoint_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> >(lhs.derived(),rhs); 00198 // } 00199 00200 // selfadjoint to dense matrix 00201 00202 namespace internal { 00203 00204 template<typename Derived1, typename Derived2, int UnrollCount, bool ClearOpposite> 00205 struct triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Upper), UnrollCount, ClearOpposite> 00206 { 00207 enum { 00208 col = (UnrollCount-1) / Derived1::RowsAtCompileTime, 00209 row = (UnrollCount-1) % Derived1::RowsAtCompileTime 00210 }; 00211 00212 static inline void run(Derived1 &dst, const Derived2 &src) 00213 { 00214 triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Upper), UnrollCount-1, ClearOpposite>::run(dst, src); 00215 00216 if(row == col) 00217 dst.coeffRef(row, col) = numext::real(src.coeff(row, col)); 00218 else if(row < col) 00219 dst.coeffRef(col, row) = numext::conj(dst.coeffRef(row, col) = src.coeff(row, col)); 00220 } 00221 }; 00222 00223 template<typename Derived1, typename Derived2, bool ClearOpposite> 00224 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Upper, 0, ClearOpposite> 00225 { 00226 static inline void run(Derived1 &, const Derived2 &) {} 00227 }; 00228 00229 template<typename Derived1, typename Derived2, int UnrollCount, bool ClearOpposite> 00230 struct triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Lower), UnrollCount, ClearOpposite> 00231 { 00232 enum { 00233 col = (UnrollCount-1) / Derived1::RowsAtCompileTime, 00234 row = (UnrollCount-1) % Derived1::RowsAtCompileTime 00235 }; 00236 00237 static inline void run(Derived1 &dst, const Derived2 &src) 00238 { 00239 triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Lower), UnrollCount-1, ClearOpposite>::run(dst, src); 00240 00241 if(row == col) 00242 dst.coeffRef(row, col) = numext::real(src.coeff(row, col)); 00243 else if(row > col) 00244 dst.coeffRef(col, row) = numext::conj(dst.coeffRef(row, col) = src.coeff(row, col)); 00245 } 00246 }; 00247 00248 template<typename Derived1, typename Derived2, bool ClearOpposite> 00249 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Lower, 0, ClearOpposite> 00250 { 00251 static inline void run(Derived1 &, const Derived2 &) {} 00252 }; 00253 00254 template<typename Derived1, typename Derived2, bool ClearOpposite> 00255 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Upper, Dynamic, ClearOpposite> 00256 { 00257 typedef typename Derived1::Index Index; 00258 static inline void run(Derived1 &dst, const Derived2 &src) 00259 { 00260 for(Index j = 0; j < dst.cols(); ++j) 00261 { 00262 for(Index i = 0; i < j; ++i) 00263 { 00264 dst.copyCoeff(i, j, src); 00265 dst.coeffRef(j,i) = numext::conj(dst.coeff(i,j)); 00266 } 00267 dst.copyCoeff(j, j, src); 00268 } 00269 } 00270 }; 00271 00272 template<typename Derived1, typename Derived2, bool ClearOpposite> 00273 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Lower, Dynamic, ClearOpposite> 00274 { 00275 static inline void run(Derived1 &dst, const Derived2 &src) 00276 { 00277 typedef typename Derived1::Index Index; 00278 for(Index i = 0; i < dst.rows(); ++i) 00279 { 00280 for(Index j = 0; j < i; ++j) 00281 { 00282 dst.copyCoeff(i, j, src); 00283 dst.coeffRef(j,i) = numext::conj(dst.coeff(i,j)); 00284 } 00285 dst.copyCoeff(i, i, src); 00286 } 00287 } 00288 }; 00289 00290 } // end namespace internal 00291 00292 /*************************************************************************** 00293 * Implementation of MatrixBase methods 00294 ***************************************************************************/ 00295 00296 template<typename Derived> 00297 template<unsigned int UpLo> 00298 typename MatrixBase<Derived>::template ConstSelfAdjointViewReturnType<UpLo>::Type 00299 MatrixBase<Derived>::selfadjointView() const 00300 { 00301 return derived(); 00302 } 00303 00304 template<typename Derived> 00305 template<unsigned int UpLo> 00306 typename MatrixBase<Derived>::template SelfAdjointViewReturnType<UpLo>::Type 00307 MatrixBase<Derived>::selfadjointView() 00308 { 00309 return derived(); 00310 } 00311 00312 } // end namespace Eigen 00313 00314 #endif // EIGEN_SELFADJOINTMATRIX_H