src/Core/SelfAdjointView.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_SELFADJOINTMATRIX_H
0011 #define EIGEN_SELFADJOINTMATRIX_H
0012
0013 namespace Eigen {
0014
0015 /** \class SelfAdjointView
0016   * \ingroup Core_Module
0017   *
0018   *
0019   * \brief Expression of a selfadjoint matrix from a triangular part of a dense matrix
0020   *
0021   * \param MatrixType the type of the dense matrix storing the coefficients
0022   * \param TriangularPart can be either \c #Lower or \c #Upper
0023   *
0024   * This class is an expression of a sefladjoint matrix from a triangular part of a matrix
0025   * with given dense storage of the coefficients. It is the return type of MatrixBase::selfadjointView()
0026   * and most of the time this is the only way that it is used.
0027   *
0028   * \sa class TriangularBase, MatrixBase::selfadjointView()
0029   */
0030
0031 namespace internal {
0032 template<typename MatrixType, unsigned int UpLo>
0033 struct traits<SelfAdjointView<MatrixType, UpLo> > : traits<MatrixType>
0034 {
0035   typedef typename ref_selector<MatrixType>::non_const_type MatrixTypeNested;
0036   typedef typename remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned;
0037   typedef MatrixType ExpressionType;
0038   typedef typename MatrixType::PlainObject FullMatrixType;
0039   enum {
0040     Mode = UpLo | SelfAdjoint,
0041     FlagsLvalueBit = is_lvalue<MatrixType>::value ? LvalueBit : 0,
0042     Flags =  MatrixTypeNestedCleaned::Flags & (HereditaryBits|FlagsLvalueBit)
0043            & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit)) // FIXME these flags should be preserved
0044   };
0045 };
0046 }
0047
0048
0049 template<typename _MatrixType, unsigned int UpLo> class SelfAdjointView
0050   : public TriangularBase<SelfAdjointView<_MatrixType, UpLo> >
0051 {
0052   public:
0053
0054     typedef _MatrixType MatrixType;
0055     typedef TriangularBase<SelfAdjointView> Base;
0056     typedef typename internal::traits<SelfAdjointView>::MatrixTypeNested MatrixTypeNested;
0057     typedef typename internal::traits<SelfAdjointView>::MatrixTypeNestedCleaned MatrixTypeNestedCleaned;
0058     typedef MatrixTypeNestedCleaned NestedExpression;
0059
0060     /** \brief The type of coefficients in this matrix */
0061     typedef typename internal::traits<SelfAdjointView>::Scalar Scalar;
0062     typedef typename MatrixType::StorageIndex StorageIndex;
0063     typedef typename internal::remove_all<typename MatrixType::ConjugateReturnType>::type MatrixConjugateReturnType;
0064     typedef SelfAdjointView<typename internal::add_const<MatrixType>::type, UpLo> ConstSelfAdjointView;
0065
0066     enum {
0067       Mode = internal::traits<SelfAdjointView>::Mode,
0068       Flags = internal::traits<SelfAdjointView>::Flags,
0069       TransposeMode = ((int(Mode) & int(Upper)) ? Lower : 0) | ((int(Mode) & int(Lower)) ? Upper : 0)
0070     };
0071     typedef typename MatrixType::PlainObject PlainObject;
0072
0073     EIGEN_DEVICE_FUNC
0074     explicit inline SelfAdjointView(MatrixType& matrix) : m_matrix(matrix)
0075     {
0076       EIGEN_STATIC_ASSERT(UpLo==Lower || UpLo==Upper,SELFADJOINTVIEW_ACCEPTS_UPPER_AND_LOWER_MODE_ONLY);
0077     }
0078
0079     EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
0080     inline Index rows() const EIGEN_NOEXCEPT { return m_matrix.rows(); }
0081     EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
0082     inline Index cols() const EIGEN_NOEXCEPT { return m_matrix.cols(); }
0083     EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
0084     inline Index outerStride() const EIGEN_NOEXCEPT { return m_matrix.outerStride(); }
0085     EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
0086     inline Index innerStride() const EIGEN_NOEXCEPT { return m_matrix.innerStride(); }
0087
0088     /** \sa MatrixBase::coeff()
0089       * \warning the coordinates must fit into the referenced triangular part
0090       */
0091     EIGEN_DEVICE_FUNC
0092     inline Scalar coeff(Index row, Index col) const
0093     {
0094       Base::check_coordinates_internal(row, col);
0095       return m_matrix.coeff(row, col);
0096     }
0097
0098     /** \sa MatrixBase::coeffRef()
0099       * \warning the coordinates must fit into the referenced triangular part
0100       */
0101     EIGEN_DEVICE_FUNC
0102     inline Scalar& coeffRef(Index row, Index col)
0103     {
0104       EIGEN_STATIC_ASSERT_LVALUE(SelfAdjointView);
0105       Base::check_coordinates_internal(row, col);
0106       return m_matrix.coeffRef(row, col);
0107     }
0108
0109     /** \internal */
0110     EIGEN_DEVICE_FUNC
0111     const MatrixTypeNestedCleaned& _expression() const { return m_matrix; }
0112
0113     EIGEN_DEVICE_FUNC
0114     const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; }
0115     EIGEN_DEVICE_FUNC
0116     MatrixTypeNestedCleaned& nestedExpression() { return m_matrix; }
0117
0118     /** Efficient triangular matrix times vector/matrix product */
0119     template<typename OtherDerived>
0120     EIGEN_DEVICE_FUNC
0121     const Product<SelfAdjointView,OtherDerived>
0122     operator*(const MatrixBase<OtherDerived>& rhs) const
0123     {
0124       return Product<SelfAdjointView,OtherDerived>(*this, rhs.derived());
0125     }
0126
0127     /** Efficient vector/matrix times triangular matrix product */
0128     template<typename OtherDerived> friend
0129     EIGEN_DEVICE_FUNC
0130     const Product<OtherDerived,SelfAdjointView>
0131     operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView& rhs)
0132     {
0133       return Product<OtherDerived,SelfAdjointView>(lhs.derived(),rhs);
0134     }
0135
0136     friend EIGEN_DEVICE_FUNC
0137     const SelfAdjointView<const EIGEN_SCALAR_BINARYOP_EXPR_RETURN_TYPE(Scalar,MatrixType,product),UpLo>
0138     operator*(const Scalar& s, const SelfAdjointView& mat)
0139     {
0140       return (s*mat.nestedExpression()).template selfadjointView<UpLo>();
0141     }
0142
0143     /** Perform a symmetric rank 2 update of the selfadjoint matrix \c *this:
0144       * \f$ this = this + \alpha u v^* + conj(\alpha) v u^* \f$
0145       * \returns a reference to \c *this
0146       *
0147       * The vectors \a u and \c v \b must be column vectors, however they can be
0148       * a adjoint expression without any overhead. Only the meaningful triangular
0149       * part of the matrix is updated, the rest is left unchanged.
0150       *
0151       * \sa rankUpdate(const MatrixBase<DerivedU>&, Scalar)
0152       */
0153     template<typename DerivedU, typename DerivedV>
0154     EIGEN_DEVICE_FUNC
0155     SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, const Scalar& alpha = Scalar(1));
0156
0157     /** Perform a symmetric rank K update of the selfadjoint matrix \c *this:
0158       * \f$ this = this + \alpha ( u u^* ) \f$ where \a u is a vector or matrix.
0159       *
0160       * \returns a reference to \c *this
0161       *
0162       * Note that to perform \f$ this = this + \alpha ( u^* u ) \f$ you can simply
0163       * call this function with u.adjoint().
0164       *
0165       * \sa rankUpdate(const MatrixBase<DerivedU>&, const MatrixBase<DerivedV>&, Scalar)
0166       */
0167     template<typename DerivedU>
0168     EIGEN_DEVICE_FUNC
0169     SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const Scalar& alpha = Scalar(1));
0170
0171     /** \returns an expression of a triangular view extracted from the current selfadjoint view of a given triangular part
0172       *
0173       * The parameter \a TriMode can have the following values: \c #Upper, \c #StrictlyUpper, \c #UnitUpper,
0174       * \c #Lower, \c #StrictlyLower, \c #UnitLower.
0175       *
0176       * If \c TriMode references the same triangular part than \c *this, then this method simply return a \c TriangularView of the nested expression,
0177       * otherwise, the nested expression is first transposed, thus returning a \c TriangularView<Transpose<MatrixType>> object.
0178       *
0179       * \sa MatrixBase::triangularView(), class TriangularView
0180       */
0181     template<unsigned int TriMode>
0182     EIGEN_DEVICE_FUNC
0183     typename internal::conditional<(TriMode&(Upper|Lower))==(UpLo&(Upper|Lower)),
0184                                    TriangularView<MatrixType,TriMode>,
0185                                    TriangularView<typename MatrixType::AdjointReturnType,TriMode> >::type
0186     triangularView() const
0187     {
0188       typename internal::conditional<(TriMode&(Upper|Lower))==(UpLo&(Upper|Lower)), MatrixType&, typename MatrixType::ConstTransposeReturnType>::type tmp1(m_matrix);
0189       typename internal::conditional<(TriMode&(Upper|Lower))==(UpLo&(Upper|Lower)), MatrixType&, typename MatrixType::AdjointReturnType>::type tmp2(tmp1);
0190       return typename internal::conditional<(TriMode&(Upper|Lower))==(UpLo&(Upper|Lower)),
0191                                    TriangularView<MatrixType,TriMode>,
0192                                    TriangularView<typename MatrixType::AdjointReturnType,TriMode> >::type(tmp2);
0193     }
0194
0195     typedef SelfAdjointView<const MatrixConjugateReturnType,UpLo> ConjugateReturnType;
0196     /** \sa MatrixBase::conjugate() const */
0197     EIGEN_DEVICE_FUNC
0198     inline const ConjugateReturnType conjugate() const
0199     { return ConjugateReturnType(m_matrix.conjugate()); }
0200
0201     /** \returns an expression of the complex conjugate of \c *this if Cond==true,
0202      *           returns \c *this otherwise.
0203      */
0204     template<bool Cond>
0205     EIGEN_DEVICE_FUNC
0206     inline typename internal::conditional<Cond,ConjugateReturnType,ConstSelfAdjointView>::type
0207     conjugateIf() const
0208     {
0209       typedef typename internal::conditional<Cond,ConjugateReturnType,ConstSelfAdjointView>::type ReturnType;
0210       return ReturnType(m_matrix.template conjugateIf<Cond>());
0211     }
0212
0213     typedef SelfAdjointView<const typename MatrixType::AdjointReturnType,TransposeMode> AdjointReturnType;
0214     /** \sa MatrixBase::adjoint() const */
0215     EIGEN_DEVICE_FUNC
0216     inline const AdjointReturnType adjoint() const
0217     { return AdjointReturnType(m_matrix.adjoint()); }
0218
0219     typedef SelfAdjointView<typename MatrixType::TransposeReturnType,TransposeMode> TransposeReturnType;
0220      /** \sa MatrixBase::transpose() */
0221     EIGEN_DEVICE_FUNC
0222     inline TransposeReturnType transpose()
0223     {
0224       EIGEN_STATIC_ASSERT_LVALUE(MatrixType)
0225       typename MatrixType::TransposeReturnType tmp(m_matrix);
0226       return TransposeReturnType(tmp);
0227     }
0228
0229     typedef SelfAdjointView<const typename MatrixType::ConstTransposeReturnType,TransposeMode> ConstTransposeReturnType;
0230     /** \sa MatrixBase::transpose() const */
0231     EIGEN_DEVICE_FUNC
0232     inline const ConstTransposeReturnType transpose() const
0233     {
0234       return ConstTransposeReturnType(m_matrix.transpose());
0235     }
0236
0237     /** \returns a const expression of the main diagonal of the matrix \c *this
0238       *
0239       * This method simply returns the diagonal of the nested expression, thus by-passing the SelfAdjointView decorator.
0240       *
0241       * \sa MatrixBase::diagonal(), class Diagonal */
0242     EIGEN_DEVICE_FUNC
0243     typename MatrixType::ConstDiagonalReturnType diagonal() const
0244     {
0245       return typename MatrixType::ConstDiagonalReturnType(m_matrix);
0246     }
0247
0248 /////////// Cholesky module ///////////
0249
0250     const LLT<PlainObject, UpLo> llt() const;
0251     const LDLT<PlainObject, UpLo> ldlt() const;
0252
0253 /////////// Eigenvalue module ///////////
0254
0255     /** Real part of #Scalar */
0256     typedef typename NumTraits<Scalar>::Real RealScalar;
0257     /** Return type of eigenvalues() */
0258     typedef Matrix<RealScalar, internal::traits<MatrixType>::ColsAtCompileTime, 1> EigenvaluesReturnType;
0259
0260     EIGEN_DEVICE_FUNC
0261     EigenvaluesReturnType eigenvalues() const;
0262     EIGEN_DEVICE_FUNC
0263     RealScalar operatorNorm() const;
0264
0265   protected:
0266     MatrixTypeNested m_matrix;
0267 };
0268
0269
0270 // template<typename OtherDerived, typename MatrixType, unsigned int UpLo>
0271 // internal::selfadjoint_matrix_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> >
0272 // operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView<MatrixType,UpLo>& rhs)
0273 // {
0274 //   return internal::matrix_selfadjoint_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> >(lhs.derived(),rhs);
0275 // }
0276
0277 // selfadjoint to dense matrix
0278
0279 namespace internal {
0280
0281 // TODO currently a selfadjoint expression has the form SelfAdjointView<.,.>
0282 //      in the future selfadjoint-ness should be defined by the expression traits
0283 //      such that Transpose<SelfAdjointView<.,.> > is valid. (currently TriangularBase::transpose() is overloaded to make it work)
0284 template<typename MatrixType, unsigned int Mode>
0285 struct evaluator_traits<SelfAdjointView<MatrixType,Mode> >
0286 {
0287   typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind;
0288   typedef SelfAdjointShape Shape;
0289 };
0290
0291 template<int UpLo, int SetOpposite, typename DstEvaluatorTypeT, typename SrcEvaluatorTypeT, typename Functor, int Version>
0292 class triangular_dense_assignment_kernel<UpLo,SelfAdjoint,SetOpposite,DstEvaluatorTypeT,SrcEvaluatorTypeT,Functor,Version>
0293   : public generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version>
0294 {
0295 protected:
0296   typedef generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version> Base;
0297   typedef typename Base::DstXprType DstXprType;
0298   typedef typename Base::SrcXprType SrcXprType;
0299   using Base::m_dst;
0300   using Base::m_src;
0301   using Base::m_functor;
0302 public:
0303
0304   typedef typename Base::DstEvaluatorType DstEvaluatorType;
0305   typedef typename Base::SrcEvaluatorType SrcEvaluatorType;
0306   typedef typename Base::Scalar Scalar;
0307   typedef typename Base::AssignmentTraits AssignmentTraits;
0308
0309
0310   EIGEN_DEVICE_FUNC triangular_dense_assignment_kernel(DstEvaluatorType &dst, const SrcEvaluatorType &src, const Functor &func, DstXprType& dstExpr)
0311     : Base(dst, src, func, dstExpr)
0312   {}
0313
0314   EIGEN_DEVICE_FUNC void assignCoeff(Index row, Index col)
0315   {
0316     eigen_internal_assert(row!=col);
0317     Scalar tmp = m_src.coeff(row,col);
0318     m_functor.assignCoeff(m_dst.coeffRef(row,col), tmp);
0319     m_functor.assignCoeff(m_dst.coeffRef(col,row), numext::conj(tmp));
0320   }
0321
0322   EIGEN_DEVICE_FUNC void assignDiagonalCoeff(Index id)
0323   {
0324     Base::assignCoeff(id,id);
0325   }
0326
0327   EIGEN_DEVICE_FUNC void assignOppositeCoeff(Index, Index)
0328   { eigen_internal_assert(false && "should never be called"); }
0329 };
0330
0331 } // end namespace internal
0332
0333 /***************************************************************************
0334 * Implementation of MatrixBase methods
0335 ***************************************************************************/
0336
0337 /** This is the const version of MatrixBase::selfadjointView() */
0338 template<typename Derived>
0339 template<unsigned int UpLo>
0340 EIGEN_DEVICE_FUNC typename MatrixBase<Derived>::template ConstSelfAdjointViewReturnType<UpLo>::Type
0341 MatrixBase<Derived>::selfadjointView() const
0342 {
0343   return typename ConstSelfAdjointViewReturnType<UpLo>::Type(derived());
0344 }
0345
0346 /** \returns an expression of a symmetric/self-adjoint view extracted from the upper or lower triangular part of the current matrix
0347   *
0348   * The parameter \a UpLo can be either \c #Upper or \c #Lower
0349   *
0350   * Example: \include MatrixBase_selfadjointView.cpp
0351   * Output: \verbinclude MatrixBase_selfadjointView.out
0352   *
0353   * \sa class SelfAdjointView
0354   */
0355 template<typename Derived>
0356 template<unsigned int UpLo>
0357 EIGEN_DEVICE_FUNC typename MatrixBase<Derived>::template SelfAdjointViewReturnType<UpLo>::Type
0358 MatrixBase<Derived>::selfadjointView()
0359 {
0360   return typename SelfAdjointViewReturnType<UpLo>::Type(derived());
0361 }
0362
0363 } // end namespace Eigen
0364
0365 #endif // EIGEN_SELFADJOINTMATRIX_H