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 namespace Eigen { 
00029 
00046 namespace internal {
00047 template<typename MatrixType, unsigned int UpLo>
00048 struct traits<SelfAdjointView<MatrixType, UpLo> > : traits<MatrixType>
00049 {
00050   typedef typename nested<MatrixType>::type MatrixTypeNested;
00051   typedef typename remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned;
00052   typedef MatrixType ExpressionType;
00053   typedef typename MatrixType::PlainObject DenseMatrixType;
00054   enum {
00055     Mode = UpLo | SelfAdjoint,
00056     Flags =  MatrixTypeNestedCleaned::Flags & (HereditaryBits)
00057            & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit)), // FIXME these flags should be preserved
00058     CoeffReadCost = MatrixTypeNestedCleaned::CoeffReadCost
00059   };
00060 };
00061 }
00062 
00063 template <typename Lhs, int LhsMode, bool LhsIsVector,
00064           typename Rhs, int RhsMode, bool RhsIsVector>
00065 struct SelfadjointProductMatrix;
00066 
00067 // FIXME could also be called SelfAdjointWrapper to be consistent with DiagonalWrapper ??
00068 template<typename MatrixType, unsigned int UpLo> class SelfAdjointView
00069   : public TriangularBase<SelfAdjointView<MatrixType, UpLo> >
00070 {
00071   public:
00072 
00073     typedef TriangularBase<SelfAdjointView> Base;
00074     typedef typename internal::traits<SelfAdjointView>::MatrixTypeNested MatrixTypeNested;
00075     typedef typename internal::traits<SelfAdjointView>::MatrixTypeNestedCleaned MatrixTypeNestedCleaned;
00076 
00078     typedef typename internal::traits<SelfAdjointView>::Scalar Scalar; 
00079 
00080     typedef typename MatrixType::Index Index;
00081 
00082     enum {
00083       Mode = internal::traits<SelfAdjointView>::Mode
00084     };
00085     typedef typename MatrixType::PlainObject PlainObject;
00086 
00087     inline SelfAdjointView(MatrixType& matrix) : m_matrix(matrix)
00088     {}
00089 
00090     inline Index rows() const { return m_matrix.rows(); }
00091     inline Index cols() const { return m_matrix.cols(); }
00092     inline Index outerStride() const { return m_matrix.outerStride(); }
00093     inline Index innerStride() const { return m_matrix.innerStride(); }
00094 
00098     inline Scalar coeff(Index row, Index col) const
00099     {
00100       Base::check_coordinates_internal(row, col);
00101       return m_matrix.coeff(row, col);
00102     }
00103 
00107     inline Scalar& coeffRef(Index row, Index col)
00108     {
00109       Base::check_coordinates_internal(row, col);
00110       return m_matrix.const_cast_derived().coeffRef(row, col);
00111     }
00112 
00114     const MatrixTypeNestedCleaned& _expression() const { return m_matrix; }
00115 
00116     const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; }
00117     MatrixTypeNestedCleaned& nestedExpression() { return *const_cast<MatrixTypeNestedCleaned*>(&m_matrix); }
00118 
00120     template<typename OtherDerived>
00121     SelfadjointProductMatrix<MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime>
00122     operator*(const MatrixBase<OtherDerived>& rhs) const
00123     {
00124       return SelfadjointProductMatrix
00125               <MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime>
00126               (m_matrix, rhs.derived());
00127     }
00128 
00130     template<typename OtherDerived> friend
00131     SelfadjointProductMatrix<OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false>
00132     operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView& rhs)
00133     {
00134       return SelfadjointProductMatrix
00135               <OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false>
00136               (lhs.derived(),rhs.m_matrix);
00137     }
00138 
00149     template<typename DerivedU, typename DerivedV>
00150     SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, Scalar alpha = Scalar(1));
00151 
00162     template<typename DerivedU>
00163     SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, Scalar alpha = Scalar(1));
00164 
00166 
00167     const LLT<PlainObject, UpLo> llt() const;
00168     const LDLT<PlainObject, UpLo> ldlt() const;
00169 
00171 
00173     typedef typename NumTraits<Scalar>::Real RealScalar;
00175     typedef Matrix<RealScalar, internal::traits<MatrixType>::ColsAtCompileTime, 1> EigenvaluesReturnType;
00176 
00177     EigenvaluesReturnType eigenvalues() const;
00178     RealScalar operatorNorm() const;
00179     
00180     #ifdef EIGEN2_SUPPORT
00181     template<typename OtherDerived>
00182     SelfAdjointView& operator=(const MatrixBase<OtherDerived>& other)
00183     {
00184       enum {
00185         OtherPart = UpLo == Upper ? StrictlyLower : StrictlyUpper
00186       };
00187       m_matrix.const_cast_derived().template triangularView<UpLo>() = other;
00188       m_matrix.const_cast_derived().template triangularView<OtherPart>() = other.adjoint();
00189       return *this;
00190     }
00191     template<typename OtherMatrixType, unsigned int OtherMode>
00192     SelfAdjointView& operator=(const TriangularView<OtherMatrixType, OtherMode>& other)
00193     {
00194       enum {
00195         OtherPart = UpLo == Upper ? StrictlyLower : StrictlyUpper
00196       };
00197       m_matrix.const_cast_derived().template triangularView<UpLo>() = other.toDenseMatrix();
00198       m_matrix.const_cast_derived().template triangularView<OtherPart>() = other.toDenseMatrix().adjoint();
00199       return *this;
00200     }
00201     #endif
00202 
00203   protected:
00204     MatrixTypeNested m_matrix;
00205 };
00206 
00207 
00208 // template<typename OtherDerived, typename MatrixType, unsigned int UpLo>
00209 // internal::selfadjoint_matrix_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> >
00210 // operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView<MatrixType,UpLo>& rhs)
00211 // {
00212 //   return internal::matrix_selfadjoint_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> >(lhs.derived(),rhs);
00213 // }
00214 
00215 // selfadjoint to dense matrix
00216 
00217 namespace internal {
00218 
00219 template<typename Derived1, typename Derived2, int UnrollCount, bool ClearOpposite>
00220 struct triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Upper), UnrollCount, ClearOpposite>
00221 {
00222   enum {
00223     col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
00224     row = (UnrollCount-1) % Derived1::RowsAtCompileTime
00225   };
00226 
00227   static inline void run(Derived1 &dst, const Derived2 &src)
00228   {
00229     triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Upper), UnrollCount-1, ClearOpposite>::run(dst, src);
00230 
00231     if(row == col)
00232       dst.coeffRef(row, col) = real(src.coeff(row, col));
00233     else if(row < col)
00234       dst.coeffRef(col, row) = conj(dst.coeffRef(row, col) = src.coeff(row, col));
00235   }
00236 };
00237 
00238 template<typename Derived1, typename Derived2, bool ClearOpposite>
00239 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Upper, 0, ClearOpposite>
00240 {
00241   static inline void run(Derived1 &, const Derived2 &) {}
00242 };
00243 
00244 template<typename Derived1, typename Derived2, int UnrollCount, bool ClearOpposite>
00245 struct triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Lower), UnrollCount, ClearOpposite>
00246 {
00247   enum {
00248     col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
00249     row = (UnrollCount-1) % Derived1::RowsAtCompileTime
00250   };
00251 
00252   static inline void run(Derived1 &dst, const Derived2 &src)
00253   {
00254     triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Lower), UnrollCount-1, ClearOpposite>::run(dst, src);
00255 
00256     if(row == col)
00257       dst.coeffRef(row, col) = real(src.coeff(row, col));
00258     else if(row > col)
00259       dst.coeffRef(col, row) = conj(dst.coeffRef(row, col) = src.coeff(row, col));
00260   }
00261 };
00262 
00263 template<typename Derived1, typename Derived2, bool ClearOpposite>
00264 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Lower, 0, ClearOpposite>
00265 {
00266   static inline void run(Derived1 &, const Derived2 &) {}
00267 };
00268 
00269 template<typename Derived1, typename Derived2, bool ClearOpposite>
00270 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Upper, Dynamic, ClearOpposite>
00271 {
00272   typedef typename Derived1::Index Index;
00273   static inline void run(Derived1 &dst, const Derived2 &src)
00274   {
00275     for(Index j = 0; j < dst.cols(); ++j)
00276     {
00277       for(Index i = 0; i < j; ++i)
00278       {
00279         dst.copyCoeff(i, j, src);
00280         dst.coeffRef(j,i) = conj(dst.coeff(i,j));
00281       }
00282       dst.copyCoeff(j, j, src);
00283     }
00284   }
00285 };
00286 
00287 template<typename Derived1, typename Derived2, bool ClearOpposite>
00288 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Lower, Dynamic, ClearOpposite>
00289 {
00290   static inline void run(Derived1 &dst, const Derived2 &src)
00291   {
00292   typedef typename Derived1::Index Index;
00293     for(Index i = 0; i < dst.rows(); ++i)
00294     {
00295       for(Index j = 0; j < i; ++j)
00296       {
00297         dst.copyCoeff(i, j, src);
00298         dst.coeffRef(j,i) = conj(dst.coeff(i,j));
00299       }
00300       dst.copyCoeff(i, i, src);
00301     }
00302   }
00303 };
00304 
00305 } // end namespace internal
00306 
00307 /***************************************************************************
00308 * Implementation of MatrixBase methods
00309 ***************************************************************************/
00310 
00311 template<typename Derived>
00312 template<unsigned int UpLo>
00313 typename MatrixBase<Derived>::template ConstSelfAdjointViewReturnType<UpLo>::Type
00314 MatrixBase<Derived>::selfadjointView() const
00315 {
00316   return derived();
00317 }
00318 
00319 template<typename Derived>
00320 template<unsigned int UpLo>
00321 typename MatrixBase<Derived>::template SelfAdjointViewReturnType<UpLo>::Type
00322 MatrixBase<Derived>::selfadjointView()
00323 {
00324   return derived();
00325 }
00326 
00327 } // end namespace Eigen
00328 
00329 #endif // EIGEN_SELFADJOINTMATRIX_H