diff options
author | Benoit Jacob <jacob.benoit.1@gmail.com> | 2010-03-08 12:37:04 -0500 |
---|---|---|
committer | Benoit Jacob <jacob.benoit.1@gmail.com> | 2010-03-08 12:37:04 -0500 |
commit | 89343a38afab4c2a677614eb928881040b763941 (patch) | |
tree | 60b9cb66d2a0491870c5681c84e5741d773d6605 | |
parent | 4293a4d1af0221ae4a9b19e28d118951ff56400b (diff) |
* Fix #97 : Householder operations on 1x1 matrices
* Fix VectorBlock on 1x1 "vectors"
* remove useless makeTrivialHouseholder function
-rw-r--r-- | Eigen/src/Core/VectorBlock.h | 21 | ||||
-rw-r--r-- | Eigen/src/Householder/Householder.h | 51 |
2 files changed, 45 insertions, 27 deletions
diff --git a/Eigen/src/Core/VectorBlock.h b/Eigen/src/Core/VectorBlock.h index 5bb7fd35d..94b56f321 100644 --- a/Eigen/src/Core/VectorBlock.h +++ b/Eigen/src/Core/VectorBlock.h @@ -59,20 +59,33 @@ template<typename VectorType, int Size> struct ei_traits<VectorBlock<VectorType, Size> > : public ei_traits<Block<VectorType, ei_traits<VectorType>::RowsAtCompileTime==1 ? 1 : Size, - ei_traits<VectorType>::ColsAtCompileTime==1 ? 1 : Size> > + ei_traits<VectorType>::ColsAtCompileTime==1 + // handle the 1x1 case. Taking a dynamic-sized vectorblock in a 1x1 xpr. + // example: when doing HouseholderQR<Matrix<float,1,1> >. + && ei_traits<VectorType>::RowsAtCompileTime!=1 + ? 1 : Size> > { }; template<typename VectorType, int Size> class VectorBlock : public Block<VectorType, ei_traits<VectorType>::RowsAtCompileTime==1 ? 1 : Size, - ei_traits<VectorType>::ColsAtCompileTime==1 ? 1 : Size> + ei_traits<VectorType>::ColsAtCompileTime==1 + // handle the 1x1 case. Taking a dynamic-sized vectorblock in a 1x1 xpr. + // example: when doing HouseholderQR<Matrix<float,1,1> >. + && ei_traits<VectorType>::RowsAtCompileTime!=1 + ? 1 : Size> { typedef Block<VectorType, ei_traits<VectorType>::RowsAtCompileTime==1 ? 1 : Size, - ei_traits<VectorType>::ColsAtCompileTime==1 ? 1 : Size> Base; + ei_traits<VectorType>::ColsAtCompileTime==1 + // handle the 1x1 case. Taking a dynamic-sized vectorblock in a 1x1 xpr. + // example: when doing HouseholderQR<Matrix<float,1,1> >. + && ei_traits<VectorType>::RowsAtCompileTime!=1 + ? 1 : Size + > Base; enum { - IsColVector = ei_traits<VectorType>::ColsAtCompileTime==1 + IsColVector = ei_traits<VectorType>::ColsAtCompileTime==1 && ei_traits<VectorType>::RowsAtCompileTime!=1 }; public: EIGEN_DENSE_PUBLIC_INTERFACE(VectorBlock) diff --git a/Eigen/src/Householder/Householder.h b/Eigen/src/Householder/Householder.h index 9d1383543..3c55e9e9c 100644 --- a/Eigen/src/Householder/Householder.h +++ b/Eigen/src/Householder/Householder.h @@ -29,19 +29,10 @@ template<int n> struct ei_decrement_size { enum { - ret = (n==1 || n==Dynamic) ? n : n-1 + ret = n==Dynamic ? n : n-1 }; }; -template<typename EssentialPart> -void makeTrivialHouseholder( - EssentialPart &essential, - typename EssentialPart::RealScalar &beta) -{ - beta = typename EssentialPart::RealScalar(0); - essential.setZero(); -} - template<typename Derived> void MatrixBase<Derived>::makeHouseholderInPlace(Scalar& tau, RealScalar& beta) { @@ -73,7 +64,7 @@ void MatrixBase<Derived>::makeHouseholder( { EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart) VectorBlock<Derived, EssentialPart::SizeAtCompileTime> tail(derived(), 1, size()-1); - + RealScalar tailSqNorm = size()==1 ? 0 : tail.squaredNorm(); Scalar c0 = coeff(0); @@ -99,12 +90,19 @@ void MatrixBase<Derived>::applyHouseholderOnTheLeft( const Scalar& tau, Scalar* workspace) { - Map<Matrix<Scalar, 1, Base::ColsAtCompileTime, PlainObject::Options, 1, Base::MaxColsAtCompileTime> > tmp(workspace,cols()); - Block<Derived, EssentialPart::SizeAtCompileTime, Derived::ColsAtCompileTime> bottom(derived(), 1, 0, rows()-1, cols()); - tmp.noalias() = essential.adjoint() * bottom; - tmp += this->row(0); - this->row(0) -= tau * tmp; - bottom.noalias() -= tau * essential * tmp; + if(rows() == 1) + { + *this *= Scalar(1)-tau; + } + else + { + Map<Matrix<Scalar, 1, Base::ColsAtCompileTime, PlainObject::Options, 1, Base::MaxColsAtCompileTime> > tmp(workspace,cols()); + Block<Derived, EssentialPart::SizeAtCompileTime, Derived::ColsAtCompileTime> bottom(derived(), 1, 0, rows()-1, cols()); + tmp.noalias() = essential.adjoint() * bottom; + tmp += this->row(0); + this->row(0) -= tau * tmp; + bottom.noalias() -= tau * essential * tmp; + } } template<typename Derived> @@ -114,12 +112,19 @@ void MatrixBase<Derived>::applyHouseholderOnTheRight( const Scalar& tau, Scalar* workspace) { - Map<Matrix<Scalar, Base::RowsAtCompileTime, 1, PlainObject::Options, Base::MaxRowsAtCompileTime, 1> > tmp(workspace,rows()); - Block<Derived, Derived::RowsAtCompileTime, EssentialPart::SizeAtCompileTime> right(derived(), 0, 1, rows(), cols()-1); - tmp.noalias() = right * essential.conjugate(); - tmp += this->col(0); - this->col(0) -= tau * tmp; - right.noalias() -= tau * tmp * essential.transpose(); + if(cols() == 1) + { + *this *= Scalar(1)-tau; + } + else + { + Map<Matrix<Scalar, Base::RowsAtCompileTime, 1, PlainObject::Options, Base::MaxRowsAtCompileTime, 1> > tmp(workspace,rows()); + Block<Derived, Derived::RowsAtCompileTime, EssentialPart::SizeAtCompileTime> right(derived(), 0, 1, rows(), cols()-1); + tmp.noalias() = right * essential.conjugate(); + tmp += this->col(0); + this->col(0) -= tau * tmp; + right.noalias() -= tau * tmp * essential.transpose(); + } } #endif // EIGEN_HOUSEHOLDER_H |