SelfAdjointView.h

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 
00015 /** \class SelfAdjointView
00016   * \ingroup Core_Module
00017   *
00018   *
00019   * \brief Expression of a selfadjoint matrix from a triangular part of a dense matrix
00020   *
00021   * \param MatrixType the type of the dense matrix storing the coefficients
00022   * \param TriangularPart can be either \c #Lower or \c #Upper
00023   *
00024   * This class is an expression of a sefladjoint matrix from a triangular part of a matrix
00025   * with given dense storage of the coefficients. It is the return type of MatrixBase::selfadjointView()
00026   * and most of the time this is the only way that it is used.
00027   *
00028   * \sa class TriangularBase, MatrixBase::selfadjointView()
00029   */
00030 
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 
00062     /** \brief The type of coefficients in this matrix */
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 
00080     /** \sa MatrixBase::coeff()
00081       * \warning the coordinates must fit into the referenced triangular part
00082       */
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 
00089     /** \sa MatrixBase::coeffRef()
00090       * \warning the coordinates must fit into the referenced triangular part
00091       */
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 
00098     /** \internal */
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 
00104     /** Efficient self-adjoint matrix times vector/matrix product */
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 
00114     /** Efficient vector/matrix times self-adjoint matrix product */
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 
00124     /** Perform a symmetric rank 2 update of the selfadjoint matrix \c *this:
00125       * \f$ this = this + \alpha u v^* + conj(\alpha) v u^* \f$
00126       * \returns a reference to \c *this
00127       *
00128       * The vectors \a u and \c v \b must be column vectors, however they can be
00129       * a adjoint expression without any overhead. Only the meaningful triangular
00130       * part of the matrix is updated, the rest is left unchanged.
00131       *
00132       * \sa rankUpdate(const MatrixBase<DerivedU>&, Scalar)
00133       */
00134     template<typename DerivedU, typename DerivedV>
00135     SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, const Scalar& alpha = Scalar(1));
00136 
00137     /** Perform a symmetric rank K update of the selfadjoint matrix \c *this:
00138       * \f$ this = this + \alpha ( u u^* ) \f$ where \a u is a vector or matrix.
00139       *
00140       * \returns a reference to \c *this
00141       *
00142       * Note that to perform \f$ this = this + \alpha ( u^* u ) \f$ you can simply
00143       * call this function with u.adjoint().
00144       *
00145       * \sa rankUpdate(const MatrixBase<DerivedU>&, const MatrixBase<DerivedV>&, Scalar)
00146       */
00147     template<typename DerivedU>
00148     SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const Scalar& alpha = Scalar(1));
00149 
00150 /////////// Cholesky module ///////////
00151 
00152     const LLT<PlainObject, UpLo> llt() const;
00153     const LDLT<PlainObject, UpLo> ldlt() const;
00154 
00155 /////////// Eigenvalue module ///////////
00156 
00157     /** Real part of #Scalar */
00158     typedef typename NumTraits<Scalar>::Real RealScalar;
00159     /** Return type of eigenvalues() */
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