Transpose.h (17606B)
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 explicit EIGEN_STRONG_INLINE Transpose(MatrixType& matrix) : m_matrix(matrix) {} 65 66 EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Transpose) 67 68 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE EIGEN_CONSTEXPR 69 Index rows() const EIGEN_NOEXCEPT { return m_matrix.cols(); } 70 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE EIGEN_CONSTEXPR 71 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 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 Index innerStride() const { return derived().nestedExpression().innerStride(); } 130 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE 131 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 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 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 DenseBase<Derived>::transpose() 183 { 184 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 DenseBase<Derived>::transpose() const 196 { 197 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 static void run(const Derived&, const OtherDerived&) 447 { 448 } 449 }; 450 451 template<typename Dst, typename Src> 452 void check_for_aliasing(const Dst &dst, const Src &src) 453 { 454 if((!Dst::IsVectorAtCompileTime) && dst.rows()>1 && dst.cols()>1) 455 internal::checkTransposeAliasing_impl<Dst, Src>::run(dst, src); 456 } 457 458 } // end namespace internal 459 460 #endif // EIGEN_NO_DEBUG 461 462 } // end namespace Eigen 463 464 #endif // EIGEN_TRANSPOSE_H