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
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044 namespace internal {
00045 template<typename MatrixType, unsigned int UpLo>
00046 struct traits<SelfAdjointView<MatrixType, UpLo> > : traits<MatrixType>
00047 {
00048 typedef typename nested<MatrixType>::type MatrixTypeNested;
00049 typedef typename remove_reference<MatrixTypeNested>::type _MatrixTypeNested;
00050 typedef MatrixType ExpressionType;
00051 enum {
00052 Mode = UpLo | SelfAdjoint,
00053 Flags = _MatrixTypeNested::Flags & (HereditaryBits)
00054 & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit)),
00055 CoeffReadCost = _MatrixTypeNested::CoeffReadCost
00056 };
00057 };
00058 }
00059
00060 template <typename Lhs, int LhsMode, bool LhsIsVector,
00061 typename Rhs, int RhsMode, bool RhsIsVector>
00062 struct SelfadjointProductMatrix;
00063
00064
00065 template<typename MatrixType, unsigned int UpLo> class SelfAdjointView
00066 : public TriangularBase<SelfAdjointView<MatrixType, UpLo> >
00067 {
00068 public:
00069
00070 typedef TriangularBase<SelfAdjointView> Base;
00071
00072
00073 typedef typename internal::traits<SelfAdjointView>::Scalar Scalar;
00074
00075 typedef typename MatrixType::Index Index;
00076
00077 enum {
00078 Mode = internal::traits<SelfAdjointView>::Mode
00079 };
00080 typedef typename MatrixType::PlainObject PlainObject;
00081
00082 inline SelfAdjointView(const MatrixType& matrix) : m_matrix(matrix)
00083 {}
00084
00085 inline Index rows() const { return m_matrix.rows(); }
00086 inline Index cols() const { return m_matrix.cols(); }
00087 inline Index outerStride() const { return m_matrix.outerStride(); }
00088 inline Index innerStride() const { return m_matrix.innerStride(); }
00089
00090
00091
00092
00093 inline Scalar coeff(Index row, Index col) const
00094 {
00095 Base::check_coordinates_internal(row, col);
00096 return m_matrix.coeff(row, col);
00097 }
00098
00099
00100
00101
00102 inline Scalar& coeffRef(Index row, Index col)
00103 {
00104 Base::check_coordinates_internal(row, col);
00105 return m_matrix.const_cast_derived().coeffRef(row, col);
00106 }
00107
00108
00109 const MatrixType& _expression() const { return m_matrix; }
00110
00111 const MatrixType& nestedExpression() const { return m_matrix; }
00112 MatrixType& nestedExpression() { return const_cast<MatrixType&>(m_matrix); }
00113
00114
00115 template<typename OtherDerived>
00116 SelfadjointProductMatrix<MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime>
00117 operator*(const MatrixBase<OtherDerived>& rhs) const
00118 {
00119 return SelfadjointProductMatrix
00120 <MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime>
00121 (m_matrix, rhs.derived());
00122 }
00123
00124
00125 template<typename OtherDerived> friend
00126 SelfadjointProductMatrix<OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false>
00127 operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView& rhs)
00128 {
00129 return SelfadjointProductMatrix
00130 <OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false>
00131 (lhs.derived(),rhs.m_matrix);
00132 }
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144 template<typename DerivedU, typename DerivedV>
00145 SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, Scalar alpha = Scalar(1));
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157 template<typename DerivedU>
00158 SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, Scalar alpha = Scalar(1));
00159
00160
00161
00162 const LLT<PlainObject, UpLo> llt() const;
00163 const LDLT<PlainObject, UpLo> ldlt() const;
00164
00165
00166
00167
00168 typedef typename NumTraits<Scalar>::Real RealScalar;
00169
00170 typedef Matrix<RealScalar, internal::traits<MatrixType>::ColsAtCompileTime, 1> EigenvaluesReturnType;
00171
00172 EigenvaluesReturnType eigenvalues() const;
00173 RealScalar operatorNorm() const;
00174
00175 protected:
00176 const typename MatrixType::Nested m_matrix;
00177 };
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189 namespace internal {
00190
00191 template<typename Derived1, typename Derived2, int UnrollCount, bool ClearOpposite>
00192 struct triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Upper), UnrollCount, ClearOpposite>
00193 {
00194 enum {
00195 col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
00196 row = (UnrollCount-1) % Derived1::RowsAtCompileTime
00197 };
00198
00199 inline static void run(Derived1 &dst, const Derived2 &src)
00200 {
00201 triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Upper), UnrollCount-1, ClearOpposite>::run(dst, src);
00202
00203 if(row == col)
00204 dst.coeffRef(row, col) = real(src.coeff(row, col));
00205 else if(row < col)
00206 dst.coeffRef(col, row) = conj(dst.coeffRef(row, col) = src.coeff(row, col));
00207 }
00208 };
00209
00210 template<typename Derived1, typename Derived2, bool ClearOpposite>
00211 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Upper, 0, ClearOpposite>
00212 {
00213 inline static void run(Derived1 &, const Derived2 &) {}
00214 };
00215
00216 template<typename Derived1, typename Derived2, int UnrollCount, bool ClearOpposite>
00217 struct triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Lower), UnrollCount, ClearOpposite>
00218 {
00219 enum {
00220 col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
00221 row = (UnrollCount-1) % Derived1::RowsAtCompileTime
00222 };
00223
00224 inline static void run(Derived1 &dst, const Derived2 &src)
00225 {
00226 triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Lower), UnrollCount-1, ClearOpposite>::run(dst, src);
00227
00228 if(row == col)
00229 dst.coeffRef(row, col) = real(src.coeff(row, col));
00230 else if(row > col)
00231 dst.coeffRef(col, row) = conj(dst.coeffRef(row, col) = src.coeff(row, col));
00232 }
00233 };
00234
00235 template<typename Derived1, typename Derived2, bool ClearOpposite>
00236 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Lower, 0, ClearOpposite>
00237 {
00238 inline static void run(Derived1 &, const Derived2 &) {}
00239 };
00240
00241 template<typename Derived1, typename Derived2, bool ClearOpposite>
00242 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Upper, Dynamic, ClearOpposite>
00243 {
00244 typedef typename Derived1::Index Index;
00245 inline static void run(Derived1 &dst, const Derived2 &src)
00246 {
00247 for(Index j = 0; j < dst.cols(); ++j)
00248 {
00249 for(Index i = 0; i < j; ++i)
00250 {
00251 dst.copyCoeff(i, j, src);
00252 dst.coeffRef(j,i) = conj(dst.coeff(i,j));
00253 }
00254 dst.copyCoeff(j, j, src);
00255 }
00256 }
00257 };
00258
00259 template<typename Derived1, typename Derived2, bool ClearOpposite>
00260 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Lower, Dynamic, ClearOpposite>
00261 {
00262 inline static void run(Derived1 &dst, const Derived2 &src)
00263 {
00264 typedef typename Derived1::Index Index;
00265 for(Index i = 0; i < dst.rows(); ++i)
00266 {
00267 for(Index j = 0; j < i; ++j)
00268 {
00269 dst.copyCoeff(i, j, src);
00270 dst.coeffRef(j,i) = conj(dst.coeff(i,j));
00271 }
00272 dst.copyCoeff(i, i, src);
00273 }
00274 }
00275 };
00276
00277 }
00278
00279
00280
00281
00282
00283 template<typename Derived>
00284 template<unsigned int UpLo>
00285 typename MatrixBase<Derived>::template ConstSelfAdjointViewReturnType<UpLo>::Type
00286 MatrixBase<Derived>::selfadjointView() const
00287 {
00288 return derived();
00289 }
00290
00291 template<typename Derived>
00292 template<unsigned int UpLo>
00293 typename MatrixBase<Derived>::template SelfAdjointViewReturnType<UpLo>::Type
00294 MatrixBase<Derived>::selfadjointView()
00295 {
00296 return derived();
00297 }
00298
00299 #endif // EIGEN_SELFADJOINTMATRIX_H