* change the nesting order of adjoint_return_type to

1 - make it easier to catch conjugate expressions
 2 - make sure there is no unecessary copy (we had NestByValue<Derived> which seems to be very bad)
* update eigensolver wrt recent changes
This commit is contained in:
Gael Guennebaud 2009-07-07 15:56:13 +02:00
parent 79877a9917
commit ea23f36c78
5 changed files with 18 additions and 12 deletions

View File

@ -116,7 +116,7 @@ template<typename MatrixType, int _UpLo> class LLT
bool m_isInitialized; bool m_isInitialized;
}; };
template<typename MatrixType/*, int UpLo*/> template<typename MatrixType>
bool ei_inplace_llt_lo(MatrixType& mat) bool ei_inplace_llt_lo(MatrixType& mat)
{ {
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
@ -157,7 +157,7 @@ bool ei_inplace_llt_lo(MatrixType& mat)
return true; return true;
} }
template<typename MatrixType/*, int UpLo*/> template<typename MatrixType>
bool ei_inplace_llt_up(MatrixType& mat) bool ei_inplace_llt_up(MatrixType& mat)
{ {
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;

View File

@ -258,8 +258,10 @@ template<typename Derived> class MatrixBase
/** \internal the return type of MatrixBase::imag() */ /** \internal the return type of MatrixBase::imag() */
typedef CwiseUnaryView<ei_scalar_imag_op<Scalar>, Derived> NonConstImagReturnType; typedef CwiseUnaryView<ei_scalar_imag_op<Scalar>, Derived> NonConstImagReturnType;
/** \internal the return type of MatrixBase::adjoint() */ /** \internal the return type of MatrixBase::adjoint() */
typedef Eigen::Transpose<NestByValue<typename ei_cleantype<ConjugateReturnType>::type> > typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex,
AdjointReturnType; CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, NestByValue<Eigen::Transpose<Derived> > >,
Transpose<Derived>
>::ret AdjointReturnType;
/** \internal the return type of MatrixBase::eigenvalues() */ /** \internal the return type of MatrixBase::eigenvalues() */
typedef Matrix<typename NumTraits<typename ei_traits<Derived>::Scalar>::Real, ei_traits<Derived>::ColsAtCompileTime, 1> EigenvaluesReturnType; typedef Matrix<typename NumTraits<typename ei_traits<Derived>::Scalar>::Real, ei_traits<Derived>::ColsAtCompileTime, 1> EigenvaluesReturnType;
/** \internal expression tyepe of a column */ /** \internal expression tyepe of a column */

View File

@ -697,7 +697,8 @@ struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCo
{}; {};
/** \internal */ /** \internal
* Overloaded to perform an efficient C += A*B */
template<typename Derived> template<typename Derived>
template<typename Lhs,typename Rhs> template<typename Lhs,typename Rhs>
inline Derived& inline Derived&
@ -710,7 +711,8 @@ MatrixBase<Derived>::operator+=(const Flagged<Product<Lhs,Rhs,CacheFriendlyProdu
return derived(); return derived();
} }
/** \internal */ /** \internal
* Overloaded to perform an efficient C -= A*B */
template<typename Derived> template<typename Derived>
template<typename Lhs,typename Rhs> template<typename Lhs,typename Rhs>
inline Derived& inline Derived&
@ -723,6 +725,8 @@ MatrixBase<Derived>::operator-=(const Flagged<Product<Lhs,Rhs,CacheFriendlyProdu
return derived(); return derived();
} }
/** \internal
* Overloaded to perform an efficient C = A*B */
template<typename Derived> template<typename Derived>
template<typename Lhs, typename Rhs> template<typename Lhs, typename Rhs>
inline Derived& MatrixBase<Derived>::lazyAssign(const Product<Lhs,Rhs,CacheFriendlyProduct>& product) inline Derived& MatrixBase<Derived>::lazyAssign(const Product<Lhs,Rhs,CacheFriendlyProduct>& product)

View File

@ -181,7 +181,7 @@ template<typename Derived>
inline const typename MatrixBase<Derived>::AdjointReturnType inline const typename MatrixBase<Derived>::AdjointReturnType
MatrixBase<Derived>::adjoint() const MatrixBase<Derived>::adjoint() const
{ {
return conjugate().nestByValue(); return transpose().nestByValue();
} }
/*************************************************************************** /***************************************************************************

View File

@ -259,11 +259,11 @@ compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors
// compute C = inv(L) A inv(L') // compute C = inv(L) A inv(L')
MatrixType matC = matA; MatrixType matC = matA;
cholB.matrixL().solveTriangularInPlace(matC); cholB.matrixL().solveInPlace(matC);
// FIXME since we currently do not support A * inv(L'), let's do (inv(L) A')' : // FIXME since we currently do not support A * inv(L'), let's do (inv(L) A')' :
matC = matC.adjoint().eval(); matC.adjointInPlace();
cholB.matrixL().template marked<LowerTriangular>().solveTriangularInPlace(matC); cholB.matrixL().solveInPlace(matC);
matC = matC.adjoint().eval(); matC.adjointInPlace();
// this version works too: // this version works too:
// matC = matC.transpose(); // matC = matC.transpose();
// cholB.matrixL().conjugate().template marked<LowerTriangular>().solveTriangularInPlace(matC); // cholB.matrixL().conjugate().template marked<LowerTriangular>().solveTriangularInPlace(matC);
@ -277,7 +277,7 @@ compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors
if (computeEigenvectors) if (computeEigenvectors)
{ {
// transform back the eigen vectors: evecs = inv(U) * evecs // transform back the eigen vectors: evecs = inv(U) * evecs
cholB.matrixL().adjoint().template marked<UpperTriangular>().solveTriangularInPlace(m_eivec); cholB.matrixU().solveInPlace(m_eivec);
for (int i=0; i<m_eivec.cols(); ++i) for (int i=0; i<m_eivec.cols(); ++i)
m_eivec.col(i) = m_eivec.col(i).normalized(); m_eivec.col(i) = m_eivec.col(i).normalized();
} }