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 // 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 
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 
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 
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 
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 
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 
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 
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 
00134     template<typename DerivedU, typename DerivedV>
00135     SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, const Scalar& alpha = Scalar(1));
00136 
00147     template<typename DerivedU>
00148     SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const Scalar& alpha = Scalar(1));
00149 
00151 
00152     const LLT<PlainObject, UpLo> llt() const;
00153     const LDLT<PlainObject, UpLo> ldlt() const;
00154 
00156 
00158     typedef typename NumTraits<Scalar>::Real RealScalar;
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


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:38:55