Branch data Line data Source code
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-2014 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 : : namespace internal {
17 : : template<typename MatrixType>
18 : : struct traits<Transpose<MatrixType> > : public traits<MatrixType>
19 : : {
20 : : typedef typename ref_selector<MatrixType>::type MatrixTypeNested;
21 : : typedef typename remove_reference<MatrixTypeNested>::type MatrixTypeNestedPlain;
22 : : enum {
23 : : RowsAtCompileTime = MatrixType::ColsAtCompileTime,
24 : : ColsAtCompileTime = MatrixType::RowsAtCompileTime,
25 : : MaxRowsAtCompileTime = MatrixType::MaxColsAtCompileTime,
26 : : MaxColsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
27 : : FlagsLvalueBit = is_lvalue<MatrixType>::value ? LvalueBit : 0,
28 : : Flags0 = traits<MatrixTypeNestedPlain>::Flags & ~(LvalueBit | NestByRefBit),
29 : : Flags1 = Flags0 | FlagsLvalueBit,
30 : : Flags = Flags1 ^ RowMajorBit,
31 : : InnerStrideAtCompileTime = inner_stride_at_compile_time<MatrixType>::ret,
32 : : OuterStrideAtCompileTime = outer_stride_at_compile_time<MatrixType>::ret
33 : : };
34 : : };
35 : : }
36 : :
37 : : template<typename MatrixType, typename StorageKind> class TransposeImpl;
38 : :
39 : : /** \class Transpose
40 : : * \ingroup Core_Module
41 : : *
42 : : * \brief Expression of the transpose of a matrix
43 : : *
44 : : * \tparam MatrixType the type of the object of which we are taking the transpose
45 : : *
46 : : * This class represents an expression of the transpose of a matrix.
47 : : * It is the return type of MatrixBase::transpose() and MatrixBase::adjoint()
48 : : * and most of the time this is the only way it is used.
49 : : *
50 : : * \sa MatrixBase::transpose(), MatrixBase::adjoint()
51 : : */
52 : : template<typename MatrixType> class Transpose
53 : : : public TransposeImpl<MatrixType,typename internal::traits<MatrixType>::StorageKind>
54 : : {
55 : : public:
56 : :
57 : : typedef typename internal::ref_selector<MatrixType>::non_const_type MatrixTypeNested;
58 : :
59 : : typedef typename TransposeImpl<MatrixType,typename internal::traits<MatrixType>::StorageKind>::Base Base;
60 : : EIGEN_GENERIC_PUBLIC_INTERFACE(Transpose)
61 : : typedef typename internal::remove_all<MatrixType>::type NestedExpression;
62 : :
63 : : EIGEN_DEVICE_FUNC
64 : 3282 : explicit EIGEN_STRONG_INLINE Transpose(MatrixType& matrix) : m_matrix(matrix) {}
65 : :
66 : 1824 : EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Transpose)
67 : :
68 : : EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE EIGEN_CONSTEXPR
69 : 8084 : Index rows() const EIGEN_NOEXCEPT { return m_matrix.cols(); }
70 : : EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE EIGEN_CONSTEXPR
71 : 9300 : Index cols() const EIGEN_NOEXCEPT { return m_matrix.rows(); }
72 : :
73 : : /** \returns the nested expression */
74 : : EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
75 : : const typename internal::remove_all<MatrixTypeNested>::type&
76 : 10122 : nestedExpression() const { return m_matrix; }
77 : :
78 : : /** \returns the nested expression */
79 : : EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
80 : : typename internal::remove_reference<MatrixTypeNested>::type&
81 : : nestedExpression() { return m_matrix; }
82 : :
83 : : /** \internal */
84 : : EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
85 : : void resize(Index nrows, Index ncols) {
86 : : m_matrix.resize(ncols,nrows);
87 : : }
88 : :
89 : : protected:
90 : : typename internal::ref_selector<MatrixType>::non_const_type m_matrix;
91 : : };
92 : :
93 : : namespace internal {
94 : :
95 : : template<typename MatrixType, bool HasDirectAccess = has_direct_access<MatrixType>::ret>
96 : : struct TransposeImpl_base
97 : : {
98 : : typedef typename dense_xpr_base<Transpose<MatrixType> >::type type;
99 : : };
100 : :
101 : : template<typename MatrixType>
102 : : struct TransposeImpl_base<MatrixType, false>
103 : : {
104 : : typedef typename dense_xpr_base<Transpose<MatrixType> >::type type;
105 : : };
106 : :
107 : : } // end namespace internal
108 : :
109 : : // Generic API dispatcher
110 : : template<typename XprType, typename StorageKind>
111 : : class TransposeImpl
112 : : : public internal::generic_xpr_base<Transpose<XprType> >::type
113 : : {
114 : : public:
115 : : typedef typename internal::generic_xpr_base<Transpose<XprType> >::type Base;
116 : : };
117 : :
118 : : template<typename MatrixType> class TransposeImpl<MatrixType,Dense>
119 : : : public internal::TransposeImpl_base<MatrixType>::type
120 : : {
121 : : public:
122 : :
123 : : typedef typename internal::TransposeImpl_base<MatrixType>::type Base;
124 : : using Base::coeffRef;
125 : : EIGEN_DENSE_PUBLIC_INTERFACE(Transpose<MatrixType>)
126 : : EIGEN_INHERIT_ASSIGNMENT_OPERATORS(TransposeImpl)
127 : :
128 : : EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
129 : 2736 : Index innerStride() const { return derived().nestedExpression().innerStride(); }
130 : : EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
131 : 2736 : Index outerStride() const { return derived().nestedExpression().outerStride(); }
132 : :
133 : : typedef typename internal::conditional<
134 : : internal::is_lvalue<MatrixType>::value,
135 : : Scalar,
136 : : const Scalar
137 : : >::type ScalarWithConstIfNotLvalue;
138 : :
139 : : EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
140 : : ScalarWithConstIfNotLvalue* data() { return derived().nestedExpression().data(); }
141 : : EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
142 : 1368 : const Scalar* data() const { return derived().nestedExpression().data(); }
143 : :
144 : : // FIXME: shall we keep the const version of coeffRef?
145 : : EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
146 : : const Scalar& coeffRef(Index rowId, Index colId) const
147 : : {
148 : : return derived().nestedExpression().coeffRef(colId, rowId);
149 : : }
150 : :
151 : : EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
152 : : const Scalar& coeffRef(Index index) const
153 : : {
154 : : return derived().nestedExpression().coeffRef(index);
155 : : }
156 : : protected:
157 : 3282 : EIGEN_DEFAULT_EMPTY_CONSTRUCTOR_AND_DESTRUCTOR(TransposeImpl)
158 : : };
159 : :
160 : : /** \returns an expression of the transpose of *this.
161 : : *
162 : : * Example: \include MatrixBase_transpose.cpp
163 : : * Output: \verbinclude MatrixBase_transpose.out
164 : : *
165 : : * \warning If you want to replace a matrix by its own transpose, do \b NOT do this:
166 : : * \code
167 : : * m = m.transpose(); // bug!!! caused by aliasing effect
168 : : * \endcode
169 : : * Instead, use the transposeInPlace() method:
170 : : * \code
171 : : * m.transposeInPlace();
172 : : * \endcode
173 : : * which gives Eigen good opportunities for optimization, or alternatively you can also do:
174 : : * \code
175 : : * m = m.transpose().eval();
176 : : * \endcode
177 : : *
178 : : * \sa transposeInPlace(), adjoint() */
179 : : template<typename Derived>
180 : : EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
181 : : Transpose<Derived>
182 : 152 : DenseBase<Derived>::transpose()
183 : : {
184 : 152 : return TransposeReturnType(derived());
185 : : }
186 : :
187 : : /** This is the const version of transpose().
188 : : *
189 : : * Make sure you read the warning for transpose() !
190 : : *
191 : : * \sa transposeInPlace(), adjoint() */
192 : : template<typename Derived>
193 : : EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
194 : : typename DenseBase<Derived>::ConstTransposeReturnType
195 : 3130 : DenseBase<Derived>::transpose() const
196 : : {
197 : 3130 : return ConstTransposeReturnType(derived());
198 : : }
199 : :
200 : : /** \returns an expression of the adjoint (i.e. conjugate transpose) of *this.
201 : : *
202 : : * Example: \include MatrixBase_adjoint.cpp
203 : : * Output: \verbinclude MatrixBase_adjoint.out
204 : : *
205 : : * \warning If you want to replace a matrix by its own adjoint, do \b NOT do this:
206 : : * \code
207 : : * m = m.adjoint(); // bug!!! caused by aliasing effect
208 : : * \endcode
209 : : * Instead, use the adjointInPlace() method:
210 : : * \code
211 : : * m.adjointInPlace();
212 : : * \endcode
213 : : * which gives Eigen good opportunities for optimization, or alternatively you can also do:
214 : : * \code
215 : : * m = m.adjoint().eval();
216 : : * \endcode
217 : : *
218 : : * \sa adjointInPlace(), transpose(), conjugate(), class Transpose, class internal::scalar_conjugate_op */
219 : : template<typename Derived>
220 : : EIGEN_DEVICE_FUNC inline const typename MatrixBase<Derived>::AdjointReturnType
221 : : MatrixBase<Derived>::adjoint() const
222 : : {
223 : : return AdjointReturnType(this->transpose());
224 : : }
225 : :
226 : : /***************************************************************************
227 : : * "in place" transpose implementation
228 : : ***************************************************************************/
229 : :
230 : : namespace internal {
231 : :
232 : : template<typename MatrixType,
233 : : bool IsSquare = (MatrixType::RowsAtCompileTime == MatrixType::ColsAtCompileTime) && MatrixType::RowsAtCompileTime!=Dynamic,
234 : : bool MatchPacketSize =
235 : : (int(MatrixType::RowsAtCompileTime) == int(internal::packet_traits<typename MatrixType::Scalar>::size))
236 : : && (internal::evaluator<MatrixType>::Flags&PacketAccessBit) >
237 : : struct inplace_transpose_selector;
238 : :
239 : : template<typename MatrixType>
240 : : struct inplace_transpose_selector<MatrixType,true,false> { // square matrix
241 : : static void run(MatrixType& m) {
242 : : m.matrix().template triangularView<StrictlyUpper>().swap(m.matrix().transpose().template triangularView<StrictlyUpper>());
243 : : }
244 : : };
245 : :
246 : : template<typename MatrixType>
247 : : struct inplace_transpose_selector<MatrixType,true,true> { // PacketSize x PacketSize
248 : : static void run(MatrixType& m) {
249 : : typedef typename MatrixType::Scalar Scalar;
250 : : typedef typename internal::packet_traits<typename MatrixType::Scalar>::type Packet;
251 : : const Index PacketSize = internal::packet_traits<Scalar>::size;
252 : : const Index Alignment = internal::evaluator<MatrixType>::Alignment;
253 : : PacketBlock<Packet> A;
254 : : for (Index i=0; i<PacketSize; ++i)
255 : : A.packet[i] = m.template packetByOuterInner<Alignment>(i,0);
256 : : internal::ptranspose(A);
257 : : for (Index i=0; i<PacketSize; ++i)
258 : : m.template writePacket<Alignment>(m.rowIndexByOuterInner(i,0), m.colIndexByOuterInner(i,0), A.packet[i]);
259 : : }
260 : : };
261 : :
262 : :
263 : : template <typename MatrixType, Index Alignment>
264 : : void BlockedInPlaceTranspose(MatrixType& m) {
265 : : typedef typename MatrixType::Scalar Scalar;
266 : : typedef typename internal::packet_traits<typename MatrixType::Scalar>::type Packet;
267 : : const Index PacketSize = internal::packet_traits<Scalar>::size;
268 : : eigen_assert(m.rows() == m.cols());
269 : : int row_start = 0;
270 : : for (; row_start + PacketSize <= m.rows(); row_start += PacketSize) {
271 : : for (int col_start = row_start; col_start + PacketSize <= m.cols(); col_start += PacketSize) {
272 : : PacketBlock<Packet> A;
273 : : if (row_start == col_start) {
274 : : for (Index i=0; i<PacketSize; ++i)
275 : : A.packet[i] = m.template packetByOuterInner<Alignment>(row_start + i,col_start);
276 : : internal::ptranspose(A);
277 : : for (Index i=0; i<PacketSize; ++i)
278 : : m.template writePacket<Alignment>(m.rowIndexByOuterInner(row_start + i, col_start), m.colIndexByOuterInner(row_start + i,col_start), A.packet[i]);
279 : : } else {
280 : : PacketBlock<Packet> B;
281 : : for (Index i=0; i<PacketSize; ++i) {
282 : : A.packet[i] = m.template packetByOuterInner<Alignment>(row_start + i,col_start);
283 : : B.packet[i] = m.template packetByOuterInner<Alignment>(col_start + i, row_start);
284 : : }
285 : : internal::ptranspose(A);
286 : : internal::ptranspose(B);
287 : : for (Index i=0; i<PacketSize; ++i) {
288 : : m.template writePacket<Alignment>(m.rowIndexByOuterInner(row_start + i, col_start), m.colIndexByOuterInner(row_start + i,col_start), B.packet[i]);
289 : : m.template writePacket<Alignment>(m.rowIndexByOuterInner(col_start + i, row_start), m.colIndexByOuterInner(col_start + i,row_start), A.packet[i]);
290 : : }
291 : : }
292 : : }
293 : : }
294 : : for (Index row = row_start; row < m.rows(); ++row) {
295 : : m.matrix().row(row).head(row).swap(
296 : : m.matrix().col(row).head(row).transpose());
297 : : }
298 : : }
299 : :
300 : : template<typename MatrixType,bool MatchPacketSize>
301 : : struct inplace_transpose_selector<MatrixType,false,MatchPacketSize> { // non square or dynamic matrix
302 : : static void run(MatrixType& m) {
303 : : typedef typename MatrixType::Scalar Scalar;
304 : : if (m.rows() == m.cols()) {
305 : : const Index PacketSize = internal::packet_traits<Scalar>::size;
306 : : if (!NumTraits<Scalar>::IsComplex && m.rows() >= PacketSize) {
307 : : if ((m.rows() % PacketSize) == 0)
308 : : BlockedInPlaceTranspose<MatrixType,internal::evaluator<MatrixType>::Alignment>(m);
309 : : else
310 : : BlockedInPlaceTranspose<MatrixType,Unaligned>(m);
311 : : }
312 : : else {
313 : : m.matrix().template triangularView<StrictlyUpper>().swap(m.matrix().transpose().template triangularView<StrictlyUpper>());
314 : : }
315 : : } else {
316 : : m = m.transpose().eval();
317 : : }
318 : : }
319 : : };
320 : :
321 : :
322 : : } // end namespace internal
323 : :
324 : : /** This is the "in place" version of transpose(): it replaces \c *this by its own transpose.
325 : : * Thus, doing
326 : : * \code
327 : : * m.transposeInPlace();
328 : : * \endcode
329 : : * has the same effect on m as doing
330 : : * \code
331 : : * m = m.transpose().eval();
332 : : * \endcode
333 : : * and is faster and also safer because in the latter line of code, forgetting the eval() results
334 : : * in a bug caused by \ref TopicAliasing "aliasing".
335 : : *
336 : : * Notice however that this method is only useful if you want to replace a matrix by its own transpose.
337 : : * If you just need the transpose of a matrix, use transpose().
338 : : *
339 : : * \note if the matrix is not square, then \c *this must be a resizable matrix.
340 : : * This excludes (non-square) fixed-size matrices, block-expressions and maps.
341 : : *
342 : : * \sa transpose(), adjoint(), adjointInPlace() */
343 : : template<typename Derived>
344 : : EIGEN_DEVICE_FUNC inline void DenseBase<Derived>::transposeInPlace()
345 : : {
346 : : eigen_assert((rows() == cols() || (RowsAtCompileTime == Dynamic && ColsAtCompileTime == Dynamic))
347 : : && "transposeInPlace() called on a non-square non-resizable matrix");
348 : : internal::inplace_transpose_selector<Derived>::run(derived());
349 : : }
350 : :
351 : : /***************************************************************************
352 : : * "in place" adjoint implementation
353 : : ***************************************************************************/
354 : :
355 : : /** This is the "in place" version of adjoint(): it replaces \c *this by its own transpose.
356 : : * Thus, doing
357 : : * \code
358 : : * m.adjointInPlace();
359 : : * \endcode
360 : : * has the same effect on m as doing
361 : : * \code
362 : : * m = m.adjoint().eval();
363 : : * \endcode
364 : : * and is faster and also safer because in the latter line of code, forgetting the eval() results
365 : : * in a bug caused by aliasing.
366 : : *
367 : : * Notice however that this method is only useful if you want to replace a matrix by its own adjoint.
368 : : * If you just need the adjoint of a matrix, use adjoint().
369 : : *
370 : : * \note if the matrix is not square, then \c *this must be a resizable matrix.
371 : : * This excludes (non-square) fixed-size matrices, block-expressions and maps.
372 : : *
373 : : * \sa transpose(), adjoint(), transposeInPlace() */
374 : : template<typename Derived>
375 : : EIGEN_DEVICE_FUNC inline void MatrixBase<Derived>::adjointInPlace()
376 : : {
377 : : derived() = adjoint().eval();
378 : : }
379 : :
380 : : #ifndef EIGEN_NO_DEBUG
381 : :
382 : : // The following is to detect aliasing problems in most common cases.
383 : :
384 : : namespace internal {
385 : :
386 : : template<bool DestIsTransposed, typename OtherDerived>
387 : : struct check_transpose_aliasing_compile_time_selector
388 : : {
389 : : enum { ret = bool(blas_traits<OtherDerived>::IsTransposed) != DestIsTransposed };
390 : : };
391 : :
392 : : template<bool DestIsTransposed, typename BinOp, typename DerivedA, typename DerivedB>
393 : : struct check_transpose_aliasing_compile_time_selector<DestIsTransposed,CwiseBinaryOp<BinOp,DerivedA,DerivedB> >
394 : : {
395 : : enum { ret = bool(blas_traits<DerivedA>::IsTransposed) != DestIsTransposed
396 : : || bool(blas_traits<DerivedB>::IsTransposed) != DestIsTransposed
397 : : };
398 : : };
399 : :
400 : : template<typename Scalar, bool DestIsTransposed, typename OtherDerived>
401 : : struct check_transpose_aliasing_run_time_selector
402 : : {
403 : : static bool run(const Scalar* dest, const OtherDerived& src)
404 : : {
405 : : return (bool(blas_traits<OtherDerived>::IsTransposed) != DestIsTransposed) && (dest!=0 && dest==(const Scalar*)extract_data(src));
406 : : }
407 : : };
408 : :
409 : : template<typename Scalar, bool DestIsTransposed, typename BinOp, typename DerivedA, typename DerivedB>
410 : : struct check_transpose_aliasing_run_time_selector<Scalar,DestIsTransposed,CwiseBinaryOp<BinOp,DerivedA,DerivedB> >
411 : : {
412 : : static bool run(const Scalar* dest, const CwiseBinaryOp<BinOp,DerivedA,DerivedB>& src)
413 : : {
414 : : return ((blas_traits<DerivedA>::IsTransposed != DestIsTransposed) && (dest!=0 && dest==(const Scalar*)extract_data(src.lhs())))
415 : : || ((blas_traits<DerivedB>::IsTransposed != DestIsTransposed) && (dest!=0 && dest==(const Scalar*)extract_data(src.rhs())));
416 : : }
417 : : };
418 : :
419 : : // the following selector, checkTransposeAliasing_impl, based on MightHaveTransposeAliasing,
420 : : // is because when the condition controlling the assert is known at compile time, ICC emits a warning.
421 : : // This is actually a good warning: in expressions that don't have any transposing, the condition is
422 : : // known at compile time to be false, and using that, we can avoid generating the code of the assert again
423 : : // and again for all these expressions that don't need it.
424 : :
425 : : template<typename Derived, typename OtherDerived,
426 : : bool MightHaveTransposeAliasing
427 : : = check_transpose_aliasing_compile_time_selector
428 : : <blas_traits<Derived>::IsTransposed,OtherDerived>::ret
429 : : >
430 : : struct checkTransposeAliasing_impl
431 : : {
432 : : static void run(const Derived& dst, const OtherDerived& other)
433 : : {
434 : : eigen_assert((!check_transpose_aliasing_run_time_selector
435 : : <typename Derived::Scalar,blas_traits<Derived>::IsTransposed,OtherDerived>
436 : : ::run(extract_data(dst), other))
437 : : && "aliasing detected during transposition, use transposeInPlace() "
438 : : "or evaluate the rhs into a temporary using .eval()");
439 : :
440 : : }
441 : : };
442 : :
443 : : template<typename Derived, typename OtherDerived>
444 : : struct checkTransposeAliasing_impl<Derived, OtherDerived, false>
445 : : {
446 : 764 : static void run(const Derived&, const OtherDerived&)
447 : : {
448 : 764 : }
449 : : };
450 : :
451 : : template<typename Dst, typename Src>
452 : 100919 : void check_for_aliasing(const Dst &dst, const Src &src)
453 : : {
454 : 1070 : if((!Dst::IsVectorAtCompileTime) && dst.rows()>1 && dst.cols()>1)
455 : 764 : internal::checkTransposeAliasing_impl<Dst, Src>::run(dst, src);
456 : 100919 : }
457 : :
458 : : } // end namespace internal
459 : :
460 : : #endif // EIGEN_NO_DEBUG
461 : :
462 : : } // end namespace Eigen
463 : :
464 : : #endif // EIGEN_TRANSPOSE_H
|
