SelfAdjointView.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_SELFADJOINTMATRIX_H
26 #define EIGEN_SELFADJOINTMATRIX_H
27 
28 namespace Eigen {
29 
46 namespace internal {
47 template<typename MatrixType, unsigned int UpLo>
48 struct traits<SelfAdjointView<MatrixType, UpLo> > : traits<MatrixType>
49 {
50  typedef typename nested<MatrixType>::type MatrixTypeNested;
51  typedef typename remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned;
52  typedef MatrixType ExpressionType;
53  typedef typename MatrixType::PlainObject DenseMatrixType;
54  enum {
55  Mode = UpLo | SelfAdjoint,
56  Flags = MatrixTypeNestedCleaned::Flags & (HereditaryBits)
57  & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit)), // FIXME these flags should be preserved
58  CoeffReadCost = MatrixTypeNestedCleaned::CoeffReadCost
59  };
60 };
61 }
62 
63 template <typename Lhs, int LhsMode, bool LhsIsVector,
64  typename Rhs, int RhsMode, bool RhsIsVector>
65 struct SelfadjointProductMatrix;
66 
67 // FIXME could also be called SelfAdjointWrapper to be consistent with DiagonalWrapper ??
68 template<typename MatrixType, unsigned int UpLo> class SelfAdjointView
69  : public TriangularBase<SelfAdjointView<MatrixType, UpLo> >
70 {
71  public:
72 
74  typedef typename internal::traits<SelfAdjointView>::MatrixTypeNested MatrixTypeNested;
75  typedef typename internal::traits<SelfAdjointView>::MatrixTypeNestedCleaned MatrixTypeNestedCleaned;
76 
78  typedef typename internal::traits<SelfAdjointView>::Scalar Scalar;
79 
80  typedef typename MatrixType::Index Index;
81 
82  enum {
83  Mode = internal::traits<SelfAdjointView>::Mode
84  };
85  typedef typename MatrixType::PlainObject PlainObject;
86 
87  inline SelfAdjointView(MatrixType& matrix) : m_matrix(matrix)
88  {}
89 
90  inline Index rows() const { return m_matrix.rows(); }
91  inline Index cols() const { return m_matrix.cols(); }
92  inline Index outerStride() const { return m_matrix.outerStride(); }
93  inline Index innerStride() const { return m_matrix.innerStride(); }
94 
98  inline Scalar coeff(Index row, Index col) const
99  {
100  Base::check_coordinates_internal(row, col);
101  return m_matrix.coeff(row, col);
102  }
103 
107  inline Scalar& coeffRef(Index row, Index col)
108  {
109  Base::check_coordinates_internal(row, col);
110  return m_matrix.const_cast_derived().coeffRef(row, col);
111  }
112 
114  const MatrixTypeNestedCleaned& _expression() const { return m_matrix; }
115 
116  const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; }
117  MatrixTypeNestedCleaned& nestedExpression() { return *const_cast<MatrixTypeNestedCleaned*>(&m_matrix); }
118 
120  template<typename OtherDerived>
121  SelfadjointProductMatrix<MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime>
123  {
124  return SelfadjointProductMatrix
125  <MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime>
126  (m_matrix, rhs.derived());
127  }
128 
130  template<typename OtherDerived> friend
131  SelfadjointProductMatrix<OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false>
133  {
134  return SelfadjointProductMatrix
135  <OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false>
136  (lhs.derived(),rhs.m_matrix);
137  }
138 
149  template<typename DerivedU, typename DerivedV>
150  SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, Scalar alpha = Scalar(1));
151 
162  template<typename DerivedU>
163  SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, Scalar alpha = Scalar(1));
164 
166 
167  const LLT<PlainObject, UpLo> llt() const;
168  const LDLT<PlainObject, UpLo> ldlt() const;
169 
171 
176 
177  EigenvaluesReturnType eigenvalues() const;
178  RealScalar operatorNorm() const;
179 
180  #ifdef EIGEN2_SUPPORT
181  template<typename OtherDerived>
182  SelfAdjointView& operator=(const MatrixBase<OtherDerived>& other)
183  {
184  enum {
185  OtherPart = UpLo == Upper ? StrictlyLower : StrictlyUpper
186  };
187  m_matrix.const_cast_derived().template triangularView<UpLo>() = other;
188  m_matrix.const_cast_derived().template triangularView<OtherPart>() = other.adjoint();
189  return *this;
190  }
191  template<typename OtherMatrixType, unsigned int OtherMode>
192  SelfAdjointView& operator=(const TriangularView<OtherMatrixType, OtherMode>& other)
193  {
194  enum {
195  OtherPart = UpLo == Upper ? StrictlyLower : StrictlyUpper
196  };
197  m_matrix.const_cast_derived().template triangularView<UpLo>() = other.toDenseMatrix();
198  m_matrix.const_cast_derived().template triangularView<OtherPart>() = other.toDenseMatrix().adjoint();
199  return *this;
200  }
201  #endif
202 
203  protected:
205 };
206 
207 
208 // template<typename OtherDerived, typename MatrixType, unsigned int UpLo>
209 // internal::selfadjoint_matrix_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> >
210 // operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView<MatrixType,UpLo>& rhs)
211 // {
212 // return internal::matrix_selfadjoint_product_returntype<OtherDerived,SelfAdjointView<MatrixType,UpLo> >(lhs.derived(),rhs);
213 // }
214 
215 // selfadjoint to dense matrix
216 
217 namespace internal {
218 
219 template<typename Derived1, typename Derived2, int UnrollCount, bool ClearOpposite>
220 struct triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Upper), UnrollCount, ClearOpposite>
221 {
222  enum {
223  col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
224  row = (UnrollCount-1) % Derived1::RowsAtCompileTime
225  };
226 
227  static inline void run(Derived1 &dst, const Derived2 &src)
228  {
229  triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Upper), UnrollCount-1, ClearOpposite>::run(dst, src);
230 
231  if(row == col)
232  dst.coeffRef(row, col) = real(src.coeff(row, col));
233  else if(row < col)
234  dst.coeffRef(col, row) = conj(dst.coeffRef(row, col) = src.coeff(row, col));
235  }
236 };
237 
238 template<typename Derived1, typename Derived2, bool ClearOpposite>
239 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Upper, 0, ClearOpposite>
240 {
241  static inline void run(Derived1 &, const Derived2 &) {}
242 };
243 
244 template<typename Derived1, typename Derived2, int UnrollCount, bool ClearOpposite>
245 struct triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Lower), UnrollCount, ClearOpposite>
246 {
247  enum {
248  col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
249  row = (UnrollCount-1) % Derived1::RowsAtCompileTime
250  };
251 
252  static inline void run(Derived1 &dst, const Derived2 &src)
253  {
254  triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Lower), UnrollCount-1, ClearOpposite>::run(dst, src);
255 
256  if(row == col)
257  dst.coeffRef(row, col) = real(src.coeff(row, col));
258  else if(row > col)
259  dst.coeffRef(col, row) = conj(dst.coeffRef(row, col) = src.coeff(row, col));
260  }
261 };
262 
263 template<typename Derived1, typename Derived2, bool ClearOpposite>
264 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Lower, 0, ClearOpposite>
265 {
266  static inline void run(Derived1 &, const Derived2 &) {}
267 };
268 
269 template<typename Derived1, typename Derived2, bool ClearOpposite>
270 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Upper, Dynamic, ClearOpposite>
271 {
272  typedef typename Derived1::Index Index;
273  static inline void run(Derived1 &dst, const Derived2 &src)
274  {
275  for(Index j = 0; j < dst.cols(); ++j)
276  {
277  for(Index i = 0; i < j; ++i)
278  {
279  dst.copyCoeff(i, j, src);
280  dst.coeffRef(j,i) = conj(dst.coeff(i,j));
281  }
282  dst.copyCoeff(j, j, src);
283  }
284  }
285 };
286 
287 template<typename Derived1, typename Derived2, bool ClearOpposite>
288 struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Lower, Dynamic, ClearOpposite>
289 {
290  static inline void run(Derived1 &dst, const Derived2 &src)
291  {
292  typedef typename Derived1::Index Index;
293  for(Index i = 0; i < dst.rows(); ++i)
294  {
295  for(Index j = 0; j < i; ++j)
296  {
297  dst.copyCoeff(i, j, src);
298  dst.coeffRef(j,i) = conj(dst.coeff(i,j));
299  }
300  dst.copyCoeff(i, i, src);
301  }
302  }
303 };
304 
305 } // end namespace internal
306 
307 /***************************************************************************
308 * Implementation of MatrixBase methods
309 ***************************************************************************/
310 
311 template<typename Derived>
312 template<unsigned int UpLo>
313 typename MatrixBase<Derived>::template ConstSelfAdjointViewReturnType<UpLo>::Type
315 {
316  return derived();
317 }
318 
319 template<typename Derived>
320 template<unsigned int UpLo>
323 {
324  return derived();
325 }
326 
327 } // end namespace Eigen
328 
329 #endif // EIGEN_SELFADJOINTMATRIX_H