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


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:21