00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef EIGEN_SELFADJOINTMATRIX_H
00026 #define EIGEN_SELFADJOINTMATRIX_H
00027
00028 namespace Eigen {
00029
00046 namespace internal {
00047 template<typename MatrixType, unsigned int UpLo>
00048 struct traits<SelfAdjointView<MatrixType, UpLo> > : traits<MatrixType>
00049 {
00050 typedef typename nested<MatrixType>::type MatrixTypeNested;
00051 typedef typename remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned;
00052 typedef MatrixType ExpressionType;
00053 typedef typename MatrixType::PlainObject DenseMatrixType;
00054 enum {
00055 Mode = UpLo | SelfAdjoint,
00056 Flags = MatrixTypeNestedCleaned::Flags & (HereditaryBits)
00057 & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit)),
00058 CoeffReadCost = MatrixTypeNestedCleaned::CoeffReadCost
00059 };
00060 };
00061 }
00062
00063 template <typename Lhs, int LhsMode, bool LhsIsVector,
00064 typename Rhs, int RhsMode, bool RhsIsVector>
00065 struct SelfadjointProductMatrix;
00066
00067
00068 template<typename MatrixType, unsigned int UpLo> class SelfAdjointView
00069 : public TriangularBase<SelfAdjointView<MatrixType, UpLo> >
00070 {
00071 public:
00072
00073 typedef TriangularBase<SelfAdjointView> Base;
00074 typedef typename internal::traits<SelfAdjointView>::MatrixTypeNested MatrixTypeNested;
00075 typedef typename internal::traits<SelfAdjointView>::MatrixTypeNestedCleaned MatrixTypeNestedCleaned;
00076
00078 typedef typename internal::traits<SelfAdjointView>::Scalar Scalar;
00079
00080 typedef typename MatrixType::Index Index;
00081
00082 enum {
00083 Mode = internal::traits<SelfAdjointView>::Mode
00084 };
00085 typedef typename MatrixType::PlainObject PlainObject;
00086
00087 inline SelfAdjointView(MatrixType& matrix) : m_matrix(matrix)
00088 {}
00089
00090 inline Index rows() const { return m_matrix.rows(); }
00091 inline Index cols() const { return m_matrix.cols(); }
00092 inline Index outerStride() const { return m_matrix.outerStride(); }
00093 inline Index innerStride() const { return m_matrix.innerStride(); }
00094
00098 inline Scalar coeff(Index row, Index col) const
00099 {
00100 Base::check_coordinates_internal(row, col);
00101 return m_matrix.coeff(row, col);
00102 }
00103
00107 inline Scalar& coeffRef(Index row, Index col)
00108 {
00109 Base::check_coordinates_internal(row, col);
00110 return m_matrix.const_cast_derived().coeffRef(row, col);
00111 }
00112
00114 const MatrixTypeNestedCleaned& _expression() const { return m_matrix; }
00115
00116 const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; }
00117 MatrixTypeNestedCleaned& nestedExpression() { return *const_cast<MatrixTypeNestedCleaned*>(&m_matrix); }
00118
00120 template<typename OtherDerived>
00121 SelfadjointProductMatrix<MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime>
00122 operator*(const MatrixBase<OtherDerived>& rhs) const
00123 {
00124 return SelfadjointProductMatrix
00125 <MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime>
00126 (m_matrix, rhs.derived());
00127 }
00128
00130 template<typename OtherDerived> friend
00131 SelfadjointProductMatrix<OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false>
00132 operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView& rhs)
00133 {
00134 return SelfadjointProductMatrix
00135 <OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false>
00136 (lhs.derived(),rhs.m_matrix);
00137 }
00138
00149 template<typename DerivedU, typename DerivedV>
00150 SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, Scalar alpha = Scalar(1));
00151
00162 template<typename DerivedU>
00163 SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, Scalar alpha = Scalar(1));
00164
00166
00167 const LLT<PlainObject, UpLo> llt() const;
00168 const LDLT<PlainObject, UpLo> ldlt() const;
00169
00171
00173 typedef typename NumTraits<Scalar>::Real RealScalar;
00175 typedef Matrix<RealScalar, internal::traits<MatrixType>::ColsAtCompileTime, 1> EigenvaluesReturnType;
00176
00177 EigenvaluesReturnType eigenvalues() const;
00178 RealScalar operatorNorm() const;
00179
00180 #ifdef EIGEN2_SUPPORT
00181 template<typename OtherDerived>
00182 SelfAdjointView& operator=(const MatrixBase<OtherDerived>& other)
00183 {
00184 enum {
00185 OtherPart = UpLo == Upper ? StrictlyLower : StrictlyUpper
00186 };
00187 m_matrix.const_cast_derived().template triangularView<UpLo>() = other;
00188 m_matrix.const_cast_derived().template triangularView<OtherPart>() = other.adjoint();
00189 return *this;
00190 }
00191 template<typename OtherMatrixType, unsigned int OtherMode>
00192 SelfAdjointView& operator=(const TriangularView<OtherMatrixType, OtherMode>& other)
00193 {
00194 enum {
00195 OtherPart = UpLo == Upper ? StrictlyLower : StrictlyUpper
00196 };
00197 m_matrix.const_cast_derived().template triangularView<UpLo>() = other.toDenseMatrix();
00198 m_matrix.const_cast_derived().template triangularView<OtherPart>() = other.toDenseMatrix().adjoint();
00199 return *this;
00200 }
00201 #endif
00202
00203 protected:
00204 MatrixTypeNested m_matrix;
00205 };
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217 namespace internal {
00218
00219 template<typename Derived1, typename Derived2, int UnrollCount, bool ClearOpposite>
00220 struct triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Upper), UnrollCount, ClearOpposite>
00221 {
00222 enum {
00223 col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
00224 row = (UnrollCount-1) % Derived1::RowsAtCompileTime
00225 };
00226
00227 static inline void run(Derived1 &dst, const Derived2 &src)
00228 {
00229 triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Upper), UnrollCount-1, ClearOpposite>::run(dst, src);
00230
00231 if(row == col)
00232 dst.coeffRef(row, col) = real(src.coeff(row, col));
00233 else if(row < col)
00234 dst.coeffRef(col, row) = conj(dst.coeffRef(row, col) = src.coeff(row, col));
00235 }
00236 };
00237
00238 template<typename Derived1, typename Derived2, bool ClearOpposite>
00239 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Upper, 0, ClearOpposite>
00240 {
00241 static inline void run(Derived1 &, const Derived2 &) {}
00242 };
00243
00244 template<typename Derived1, typename Derived2, int UnrollCount, bool ClearOpposite>
00245 struct triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Lower), UnrollCount, ClearOpposite>
00246 {
00247 enum {
00248 col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
00249 row = (UnrollCount-1) % Derived1::RowsAtCompileTime
00250 };
00251
00252 static inline void run(Derived1 &dst, const Derived2 &src)
00253 {
00254 triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Lower), UnrollCount-1, ClearOpposite>::run(dst, src);
00255
00256 if(row == col)
00257 dst.coeffRef(row, col) = real(src.coeff(row, col));
00258 else if(row > col)
00259 dst.coeffRef(col, row) = conj(dst.coeffRef(row, col) = src.coeff(row, col));
00260 }
00261 };
00262
00263 template<typename Derived1, typename Derived2, bool ClearOpposite>
00264 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Lower, 0, ClearOpposite>
00265 {
00266 static inline void run(Derived1 &, const Derived2 &) {}
00267 };
00268
00269 template<typename Derived1, typename Derived2, bool ClearOpposite>
00270 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Upper, Dynamic, ClearOpposite>
00271 {
00272 typedef typename Derived1::Index Index;
00273 static inline void run(Derived1 &dst, const Derived2 &src)
00274 {
00275 for(Index j = 0; j < dst.cols(); ++j)
00276 {
00277 for(Index i = 0; i < j; ++i)
00278 {
00279 dst.copyCoeff(i, j, src);
00280 dst.coeffRef(j,i) = conj(dst.coeff(i,j));
00281 }
00282 dst.copyCoeff(j, j, src);
00283 }
00284 }
00285 };
00286
00287 template<typename Derived1, typename Derived2, bool ClearOpposite>
00288 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Lower, Dynamic, ClearOpposite>
00289 {
00290 static inline void run(Derived1 &dst, const Derived2 &src)
00291 {
00292 typedef typename Derived1::Index Index;
00293 for(Index i = 0; i < dst.rows(); ++i)
00294 {
00295 for(Index j = 0; j < i; ++j)
00296 {
00297 dst.copyCoeff(i, j, src);
00298 dst.coeffRef(j,i) = conj(dst.coeff(i,j));
00299 }
00300 dst.copyCoeff(i, i, src);
00301 }
00302 }
00303 };
00304
00305 }
00306
00307
00308
00309
00310
00311 template<typename Derived>
00312 template<unsigned int UpLo>
00313 typename MatrixBase<Derived>::template ConstSelfAdjointViewReturnType<UpLo>::Type
00314 MatrixBase<Derived>::selfadjointView() const
00315 {
00316 return derived();
00317 }
00318
00319 template<typename Derived>
00320 template<unsigned int UpLo>
00321 typename MatrixBase<Derived>::template SelfAdjointViewReturnType<UpLo>::Type
00322 MatrixBase<Derived>::selfadjointView()
00323 {
00324 return derived();
00325 }
00326
00327 }
00328
00329 #endif // EIGEN_SELFADJOINTMATRIX_H