1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
5 // Copyright (C) 2009-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11 #ifndef EIGEN_TRANSPOSE_H
12 #define EIGEN_TRANSPOSE_H
13 
14 namespace Eigen {
15 
16 /** \class Transpose
17   * \ingroup Core_Module
18   *
19   * \brief Expression of the transpose of a matrix
20   *
21   * \param MatrixType the type of the object of which we are taking the transpose
22   *
23   * This class represents an expression of the transpose of a matrix.
24   * It is the return type of MatrixBase::transpose() and MatrixBase::adjoint()
25   * and most of the time this is the only way it is used.
26   *
27   * \sa MatrixBase::transpose(), MatrixBase::adjoint()
28   */
29 
30 namespace internal {
31 template<typename MatrixType>
32 struct traits<Transpose<MatrixType> > : traits<MatrixType>
33 {
34   typedef typename MatrixType::Scalar Scalar;
35   typedef typename nested<MatrixType>::type MatrixTypeNested;
36   typedef typename remove_reference<MatrixTypeNested>::type MatrixTypeNestedPlain;
37   typedef typename traits<MatrixType>::StorageKind StorageKind;
38   typedef typename traits<MatrixType>::XprKind XprKind;
39   enum {
40     RowsAtCompileTime = MatrixType::ColsAtCompileTime,
41     ColsAtCompileTime = MatrixType::RowsAtCompileTime,
42     MaxRowsAtCompileTime = MatrixType::MaxColsAtCompileTime,
43     MaxColsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
44     FlagsLvalueBit = is_lvalue<MatrixType>::value ? LvalueBit : 0,
45     Flags0 = MatrixTypeNestedPlain::Flags & ~(LvalueBit | NestByRefBit),
46     Flags1 = Flags0 | FlagsLvalueBit,
47     Flags = Flags1 ^ RowMajorBit,
48     CoeffReadCost = MatrixTypeNestedPlain::CoeffReadCost,
49     InnerStrideAtCompileTime = inner_stride_at_compile_time<MatrixType>::ret,
50     OuterStrideAtCompileTime = outer_stride_at_compile_time<MatrixType>::ret
51   };
52 };
53 }
54 
55 template<typename MatrixType, typename StorageKind> class TransposeImpl;
56 
57 template<typename MatrixType> class Transpose
58   : public TransposeImpl<MatrixType,typename internal::traits<MatrixType>::StorageKind>
59 {
60   public:
61 
62     typedef typename TransposeImpl<MatrixType,typename internal::traits<MatrixType>::StorageKind>::Base Base;
63     EIGEN_GENERIC_PUBLIC_INTERFACE(Transpose)
64 
65     inline Transpose(MatrixType& a_matrix) : m_matrix(a_matrix) {}
66 
67     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Transpose)
68 
69     inline Index rows() const { return m_matrix.cols(); }
70     inline Index cols() const { return m_matrix.rows(); }
71 
72     /** \returns the nested expression */
73     const typename internal::remove_all<typename MatrixType::Nested>::type&
74     nestedExpression() const { return m_matrix; }
75 
76     /** \returns the nested expression */
77     typename internal::remove_all<typename MatrixType::Nested>::type&
78     nestedExpression() { return m_matrix.const_cast_derived(); }
79 
80   protected:
81     typename MatrixType::Nested m_matrix;
82 };
83 
84 namespace internal {
85 
86 template<typename MatrixType, bool HasDirectAccess = has_direct_access<MatrixType>::ret>
87 struct TransposeImpl_base
88 {
89   typedef typename dense_xpr_base<Transpose<MatrixType> >::type type;
90 };
91 
92 template<typename MatrixType>
93 struct TransposeImpl_base<MatrixType, false>
94 {
95   typedef typename dense_xpr_base<Transpose<MatrixType> >::type type;
96 };
97 
98 } // end namespace internal
99 
100 template<typename MatrixType> class TransposeImpl<MatrixType,Dense>
101   : public internal::TransposeImpl_base<MatrixType>::type
102 {
103   public:
104 
105     typedef typename internal::TransposeImpl_base<MatrixType>::type Base;
106     EIGEN_DENSE_PUBLIC_INTERFACE(Transpose<MatrixType>)
107     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(TransposeImpl)
108 
109     inline Index innerStride() const { return derived().nestedExpression().innerStride(); }
110     inline Index outerStride() const { return derived().nestedExpression().outerStride(); }
111 
112     typedef typename internal::conditional<
113                        internal::is_lvalue<MatrixType>::value,
114                        Scalar,
115                        const Scalar
116                      >::type ScalarWithConstIfNotLvalue;
117 
118     inline ScalarWithConstIfNotLvalue* data() { return derived().nestedExpression().data(); }
119     inline const Scalar* data() const { return derived().nestedExpression().data(); }
120 
121     inline ScalarWithConstIfNotLvalue& coeffRef(Index rowId, Index colId)
122     {
123       EIGEN_STATIC_ASSERT_LVALUE(MatrixType)
124       return derived().nestedExpression().const_cast_derived().coeffRef(colId, rowId);
125     }
126 
127     inline ScalarWithConstIfNotLvalue& coeffRef(Index index)
128     {
129       EIGEN_STATIC_ASSERT_LVALUE(MatrixType)
130       return derived().nestedExpression().const_cast_derived().coeffRef(index);
131     }
132 
133     inline const Scalar& coeffRef(Index rowId, Index colId) const
134     {
135       return derived().nestedExpression().coeffRef(colId, rowId);
136     }
137 
138     inline const Scalar& coeffRef(Index index) const
139     {
140       return derived().nestedExpression().coeffRef(index);
141     }
142 
143     inline CoeffReturnType coeff(Index rowId, Index colId) const
144     {
145       return derived().nestedExpression().coeff(colId, rowId);
146     }
147 
148     inline CoeffReturnType coeff(Index index) const
149     {
150       return derived().nestedExpression().coeff(index);
151     }
152 
153     template<int LoadMode>
154     inline const PacketScalar packet(Index rowId, Index colId) const
155     {
156       return derived().nestedExpression().template packet<LoadMode>(colId, rowId);
157     }
158 
159     template<int LoadMode>
160     inline void writePacket(Index rowId, Index colId, const PacketScalar& x)
161     {
162       derived().nestedExpression().const_cast_derived().template writePacket<LoadMode>(colId, rowId, x);
163     }
164 
165     template<int LoadMode>
166     inline const PacketScalar packet(Index index) const
167     {
168       return derived().nestedExpression().template packet<LoadMode>(index);
169     }
170 
171     template<int LoadMode>
172     inline void writePacket(Index index, const PacketScalar& x)
173     {
174       derived().nestedExpression().const_cast_derived().template writePacket<LoadMode>(index, x);
175     }
176 };
177 
178 /** \returns an expression of the transpose of *this.
179   *
180   * Example: \include MatrixBase_transpose.cpp
181   * Output: \verbinclude MatrixBase_transpose.out
182   *
183   * \warning If you want to replace a matrix by its own transpose, do \b NOT do this:
184   * \code
185   * m = m.transpose(); // bug!!! caused by aliasing effect
186   * \endcode
187   * Instead, use the transposeInPlace() method:
188   * \code
189   * m.transposeInPlace();
190   * \endcode
191   * which gives Eigen good opportunities for optimization, or alternatively you can also do:
192   * \code
193   * m = m.transpose().eval();
194   * \endcode
195   *
196   * \sa transposeInPlace(), adjoint() */
197 template<typename Derived>
198 inline Transpose<Derived>
199 DenseBase<Derived>::transpose()
200 {
201   return derived();
202 }
203 
204 /** This is the const version of transpose().
205   *
206   * Make sure you read the warning for transpose() !
207   *
208   * \sa transposeInPlace(), adjoint() */
209 template<typename Derived>
210 inline typename DenseBase<Derived>::ConstTransposeReturnType
211 DenseBase<Derived>::transpose() const
212 {
213   return ConstTransposeReturnType(derived());
214 }
215 
216 /** \returns an expression of the adjoint (i.e. conjugate transpose) of *this.
217   *
218   * Example: \include MatrixBase_adjoint.cpp
219   * Output: \verbinclude MatrixBase_adjoint.out
220   *
221   * \warning If you want to replace a matrix by its own adjoint, do \b NOT do this:
222   * \code
223   * m = m.adjoint(); // bug!!! caused by aliasing effect
224   * \endcode
225   * Instead, use the adjointInPlace() method:
226   * \code
227   * m.adjointInPlace();
228   * \endcode
229   * which gives Eigen good opportunities for optimization, or alternatively you can also do:
230   * \code
231   * m = m.adjoint().eval();
232   * \endcode
233   *
234   * \sa adjointInPlace(), transpose(), conjugate(), class Transpose, class internal::scalar_conjugate_op */
235 template<typename Derived>
236 inline const typename MatrixBase<Derived>::AdjointReturnType
237 MatrixBase<Derived>::adjoint() const
238 {
239   return this->transpose(); // in the complex case, the .conjugate() is be implicit here
240                             // due to implicit conversion to return type
241 }
242 
243 /***************************************************************************
244 * "in place" transpose implementation
245 ***************************************************************************/
246 
247 namespace internal {
248 
249 template<typename MatrixType,
250   bool IsSquare = (MatrixType::RowsAtCompileTime == MatrixType::ColsAtCompileTime) && MatrixType::RowsAtCompileTime!=Dynamic>
251 struct inplace_transpose_selector;
252 
253 template<typename MatrixType>
254 struct inplace_transpose_selector<MatrixType,true> { // square matrix
255   static void run(MatrixType& m) {
256     m.matrix().template triangularView<StrictlyUpper>().swap(m.matrix().transpose());
257   }
258 };
259 
260 template<typename MatrixType>
261 struct inplace_transpose_selector<MatrixType,false> { // non square matrix
262   static void run(MatrixType& m) {
263     if (m.rows()==m.cols())
264       m.matrix().template triangularView<StrictlyUpper>().swap(m.matrix().transpose());
265     else
266       m = m.transpose().eval();
267   }
268 };
269 
270 } // end namespace internal
271 
272 /** This is the "in place" version of transpose(): it replaces \c *this by its own transpose.
273   * Thus, doing
274   * \code
275   * m.transposeInPlace();
276   * \endcode
277   * has the same effect on m as doing
278   * \code
279   * m = m.transpose().eval();
280   * \endcode
281   * and is faster and also safer because in the latter line of code, forgetting the eval() results
282   * in a bug caused by \ref TopicAliasing "aliasing".
283   *
284   * Notice however that this method is only useful if you want to replace a matrix by its own transpose.
285   * If you just need the transpose of a matrix, use transpose().
286   *
287   * \note if the matrix is not square, then \c *this must be a resizable matrix.
288   * This excludes (non-square) fixed-size matrices, block-expressions and maps.
289   *
290   * \sa transpose(), adjoint(), adjointInPlace() */
291 template<typename Derived>
292 inline void DenseBase<Derived>::transposeInPlace()
293 {
294   eigen_assert((rows() == cols() || (RowsAtCompileTime == Dynamic && ColsAtCompileTime == Dynamic))
295                && "transposeInPlace() called on a non-square non-resizable matrix");
296   internal::inplace_transpose_selector<Derived>::run(derived());
297 }
298 
299 /***************************************************************************
300 * "in place" adjoint implementation
301 ***************************************************************************/
302 
303 /** This is the "in place" version of adjoint(): it replaces \c *this by its own transpose.
304   * Thus, doing
305   * \code
306   * m.adjointInPlace();
307   * \endcode
308   * has the same effect on m as doing
309   * \code
310   * m = m.adjoint().eval();
311   * \endcode
312   * and is faster and also safer because in the latter line of code, forgetting the eval() results
313   * in a bug caused by aliasing.
314   *
315   * Notice however that this method is only useful if you want to replace a matrix by its own adjoint.
316   * If you just need the adjoint of a matrix, use adjoint().
317   *
318   * \note if the matrix is not square, then \c *this must be a resizable matrix.
319   * This excludes (non-square) fixed-size matrices, block-expressions and maps.
320   *
321   * \sa transpose(), adjoint(), transposeInPlace() */
322 template<typename Derived>
323 inline void MatrixBase<Derived>::adjointInPlace()
324 {
325   derived() = adjoint().eval();
326 }
327 
328 #ifndef EIGEN_NO_DEBUG
329 
330 // The following is to detect aliasing problems in most common cases.
331 
332 namespace internal {
333 
334 template<typename BinOp,typename NestedXpr,typename Rhs>
335 struct blas_traits<SelfCwiseBinaryOp<BinOp,NestedXpr,Rhs> >
336  : blas_traits<NestedXpr>
337 {
338   typedef SelfCwiseBinaryOp<BinOp,NestedXpr,Rhs> XprType;
339   static inline const XprType extract(const XprType& x) { return x; }
340 };
341 
342 template<bool DestIsTransposed, typename OtherDerived>
343 struct check_transpose_aliasing_compile_time_selector
344 {
345   enum { ret = bool(blas_traits<OtherDerived>::IsTransposed) != DestIsTransposed };
346 };
347 
348 template<bool DestIsTransposed, typename BinOp, typename DerivedA, typename DerivedB>
349 struct check_transpose_aliasing_compile_time_selector<DestIsTransposed,CwiseBinaryOp<BinOp,DerivedA,DerivedB> >
350 {
351   enum { ret =    bool(blas_traits<DerivedA>::IsTransposed) != DestIsTransposed
352                || bool(blas_traits<DerivedB>::IsTransposed) != DestIsTransposed
353   };
354 };
355 
356 template<typename Scalar, bool DestIsTransposed, typename OtherDerived>
357 struct check_transpose_aliasing_run_time_selector
358 {
359   static bool run(const Scalar* dest, const OtherDerived& src)
360   {
361     return (bool(blas_traits<OtherDerived>::IsTransposed) != DestIsTransposed) && (dest!=0 && dest==(const Scalar*)extract_data(src));
362   }
363 };
364 
365 template<typename Scalar, bool DestIsTransposed, typename BinOp, typename DerivedA, typename DerivedB>
366 struct check_transpose_aliasing_run_time_selector<Scalar,DestIsTransposed,CwiseBinaryOp<BinOp,DerivedA,DerivedB> >
367 {
368   static bool run(const Scalar* dest, const CwiseBinaryOp<BinOp,DerivedA,DerivedB>& src)
369   {
370     return ((blas_traits<DerivedA>::IsTransposed != DestIsTransposed) && (dest!=0 && dest==(const Scalar*)extract_data(src.lhs())))
371         || ((blas_traits<DerivedB>::IsTransposed != DestIsTransposed) && (dest!=0 && dest==(const Scalar*)extract_data(src.rhs())));
372   }
373 };
374 
375 // the following selector, checkTransposeAliasing_impl, based on MightHaveTransposeAliasing,
376 // is because when the condition controlling the assert is known at compile time, ICC emits a warning.
377 // This is actually a good warning: in expressions that don't have any transposing, the condition is
378 // known at compile time to be false, and using that, we can avoid generating the code of the assert again
379 // and again for all these expressions that don't need it.
380 
381 template<typename Derived, typename OtherDerived,
382          bool MightHaveTransposeAliasing
383                  = check_transpose_aliasing_compile_time_selector
384                      <blas_traits<Derived>::IsTransposed,OtherDerived>::ret
385         >
386 struct checkTransposeAliasing_impl
387 {
388     static void run(const Derived& dst, const OtherDerived& other)
389     {
390         eigen_assert((!check_transpose_aliasing_run_time_selector
391                       <typename Derived::Scalar,blas_traits<Derived>::IsTransposed,OtherDerived>
392                       ::run(extract_data(dst), other))
393           && "aliasing detected during transposition, use transposeInPlace() "
394              "or evaluate the rhs into a temporary using .eval()");
395 
396     }
397 };
398 
399 template<typename Derived, typename OtherDerived>
400 struct checkTransposeAliasing_impl<Derived, OtherDerived, false>
401 {
402     static void run(const Derived&, const OtherDerived&)
403     {
404     }
405 };
406 
407 } // end namespace internal
408 
409 template<typename Derived>
410 template<typename OtherDerived>
411 void DenseBase<Derived>::checkTransposeAliasing(const OtherDerived& other) const
412 {
413     internal::checkTransposeAliasing_impl<Derived, OtherDerived>::run(derived(), other);
414 }
415 #endif
416 
417 } // end namespace Eigen
418 
419 #endif // EIGEN_TRANSPOSE_H
420