GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-14 20:56:00 +08:00
Code Issues Packages Projects Releases Wiki Activity

* Fix #97 : Householder operations on 1x1 matrices

Browse Source 
* Fix VectorBlock on 1x1 "vectors"
* remove useless makeTrivialHouseholder function
This commit is contained in:
Benoit Jacob 2010-03-08 12:37:04 -05:00
parent 4293a4d1af
commit 89343a38af
2 changed files with 45 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
Eigen/src/Core/VectorBlock.h
View File

@ -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)

51
Eigen/src/Householder/Householder.h
View File

@ -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