Main MRPT website > C++ reference
MRPT logo

SelfAdjointView.h

Go to the documentation of this file.
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 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_SELFADJOINTMATRIX_H
00026 #define EIGEN_SELFADJOINTMATRIX_H
00027 
00028 /** \class SelfAdjointView
00029   * \ingroup Core_Module
00030   *
00031   *
00032   * \brief Expression of a selfadjoint matrix from a triangular part of a dense matrix
00033   *
00034   * \param MatrixType the type of the dense matrix storing the coefficients
00035   * \param TriangularPart can be either \c Lower or \c Upper
00036   *
00037   * This class is an expression of a sefladjoint matrix from a triangular part of a matrix
00038   * with given dense storage of the coefficients. It is the return type of MatrixBase::selfadjointView()
00039   * and most of the time this is the only way that it is used.
00040   *
00041   * \sa class TriangularBase, MatrixBase::selfAdjointView()
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)), // FIXME these flags should be preserved
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 // FIXME could also be called SelfAdjointWrapper to be consistent with DiagonalWrapper ??
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     /** \brief The type of coefficients in this matrix */
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     /** \sa MatrixBase::coeff()
00091       * \warning the coordinates must fit into the referenced triangular part
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     /** \sa MatrixBase::coeffRef()
00100       * \warning the coordinates must fit into the referenced triangular part
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     /** \internal */
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     /** Efficient self-adjoint matrix times vector/matrix product */
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     /** Efficient vector/matrix times self-adjoint matrix product */
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     /** Perform a symmetric rank 2 update of the selfadjoint matrix \c *this:
00135       * \f$ this = this + \alpha u v^* + conj(\alpha) v u^* \f$
00136       * \returns a reference to \c *this
00137       *
00138       * The vectors \a u and \c v \b must be column vectors, however they can be
00139       * a adjoint expression without any overhead. Only the meaningful triangular
00140       * part of the matrix is updated, the rest is left unchanged.
00141       *
00142       * \sa rankUpdate(const MatrixBase<DerivedU>&, Scalar)
00143       */
00144     template<typename DerivedU, typename DerivedV>
00145     SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, Scalar alpha = Scalar(1));
00146 
00147     /** Perform a symmetric rank K update of the selfadjoint matrix \c *this:
00148       * \f$ this = this + \alpha ( u u^* ) \f$ where \a u is a vector or matrix.
00149       *
00150       * \returns a reference to \c *this
00151       *
00152       * Note that to perform \f$ this = this + \alpha ( u^* u ) \f$ you can simply
00153       * call this function with u.adjoint().
00154       *
00155       * \sa rankUpdate(const MatrixBase<DerivedU>&, const MatrixBase<DerivedV>&, Scalar)
00156       */
00157     template<typename DerivedU>
00158     SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, Scalar alpha = Scalar(1));
00159 
00160 /////////// Cholesky module ///////////
00161 
00162     const LLT<PlainObject, UpLo> llt() const;
00163     const LDLT<PlainObject, UpLo> ldlt() const;
00164 
00165 /////////// Eigenvalue module ///////////
00166 
00167     /** Real part of #Scalar */
00168     typedef typename NumTraits<Scalar>::Real RealScalar;
00169     /** Return type of eigenvalues() */
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 // template<typename OtherDerived, typename MatrixType, unsigned int UpLo>
00181 // internal::selfadjoint_matrix_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> >
00182 // operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView<MatrixType,UpLo>& rhs)
00183 // {
00184 //   return internal::matrix_selfadjoint_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> >(lhs.derived(),rhs);
00185 // }
00186 
00187 // selfadjoint to dense matrix
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 } // end namespace internal
00278 
00279 /***************************************************************************
00280 * Implementation of MatrixBase methods
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



Page generated by Doxygen 1.7.3 for MRPT 0.9.4 SVN: at Sat Mar 26 06:40:17 UTC 2011