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

*adapt Householder to the convention that we now favor refs over ptrs for output. Keep "workspace" as pointer because it is an array (which is now more obvious).

Browse Source 
*rename makeHouseholderSequence to householderSequence, because that's what it returns.
This commit is contained in:
Benoit Jacob 2009-11-10 21:22:20 -05:00
parent e06aa749a4
commit a440385b41
8 changed files with 28 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
Eigen/src/Core/MatrixBase.h
View File

@ -783,10 +783,10 @@ template<typename Derived> class MatrixBase
////////// Householder module /////////// ////////// Householder module ///////////
void makeHouseholderInPlace(Scalar *tau, RealScalar *beta); void makeHouseholderInPlace(Scalar& tau, RealScalar& beta);
template<typename EssentialPart> template<typename EssentialPart>
void makeHouseholder(EssentialPart *essential, void makeHouseholder(EssentialPart& essential,
Scalar *tau, RealScalar *beta) const; Scalar& tau, RealScalar& beta) const;
template<typename EssentialPart> template<typename EssentialPart>
void applyHouseholderOnTheLeft(const EssentialPart& essential, void applyHouseholderOnTheLeft(const EssentialPart& essential,
const Scalar& tau, const Scalar& tau,

2
Eigen/src/Eigenvalues/Tridiagonalization.h
View File

@ -204,7 +204,7 @@ void Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType&
int remainingSize = n-i-1; int remainingSize = n-i-1;
RealScalar beta; RealScalar beta;
Scalar h; Scalar h;
matA.col(i).end(remainingSize).makeHouseholderInPlace(&h, &beta); matA.col(i).end(remainingSize).makeHouseholderInPlace(h, beta);
// Apply similarity transformation to remaining columns, // Apply similarity transformation to remaining columns,
// i.e., A = H A H' where H = I - h v v' and v = matA.col(i).end(n-i-1) // i.e., A = H A H' where H = I - h v v' and v = matA.col(i).end(n-i-1)

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

@ -35,18 +35,18 @@ template<int n> struct ei_decrement_size
template<typename EssentialPart> template<typename EssentialPart>
void makeTrivialHouseholder( void makeTrivialHouseholder(
EssentialPart *essential, EssentialPart &essential,
typename EssentialPart::RealScalar *beta) typename EssentialPart::RealScalar &beta)
{ {
*beta = typename EssentialPart::RealScalar(0); beta = typename EssentialPart::RealScalar(0);
essential->setZero(); essential.setZero();
} }
template<typename Derived> template<typename Derived>
void MatrixBase<Derived>::makeHouseholderInPlace(Scalar *tau, RealScalar *beta) void MatrixBase<Derived>::makeHouseholderInPlace(Scalar& tau, RealScalar& beta)
{ {
VectorBlock<Derived, ei_decrement_size<SizeAtCompileTime>::ret> essentialPart(derived(), 1, size()-1); VectorBlock<Derived, ei_decrement_size<SizeAtCompileTime>::ret> essentialPart(derived(), 1, size()-1);
makeHouseholder(&essentialPart, tau, beta); makeHouseholder(essentialPart, tau, beta);
} }
/** Computes the elementary reflector H such that: /** Computes the elementary reflector H such that:
@ -67,9 +67,9 @@ void MatrixBase<Derived>::makeHouseholderInPlace(Scalar *tau, RealScalar *beta)
template<typename Derived> template<typename Derived>
template<typename EssentialPart> template<typename EssentialPart>
void MatrixBase<Derived>::makeHouseholder( void MatrixBase<Derived>::makeHouseholder(
EssentialPart *essential, EssentialPart& essential,
Scalar *tau, Scalar& tau,
RealScalar *beta) const RealScalar& beta) const
{ {
EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart) EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart)
VectorBlock<Derived, EssentialPart::SizeAtCompileTime> tail(derived(), 1, size()-1); VectorBlock<Derived, EssentialPart::SizeAtCompileTime> tail(derived(), 1, size()-1);
@ -79,16 +79,16 @@ void MatrixBase<Derived>::makeHouseholder(
if(tailSqNorm == RealScalar(0) && ei_imag(c0)==RealScalar(0)) if(tailSqNorm == RealScalar(0) && ei_imag(c0)==RealScalar(0))
{ {
*tau = 0; tau = 0;
*beta = ei_real(c0); beta = ei_real(c0);
} }
else else
{ {
*beta = ei_sqrt(ei_abs2(c0) + tailSqNorm); beta = ei_sqrt(ei_abs2(c0) + tailSqNorm);
if (ei_real(c0)>=0.) if (ei_real(c0)>=0.)
*beta = -*beta; beta = -beta;
*essential = tail / (c0 - *beta); essential = tail / (c0 - beta);
*tau = ei_conj((*beta - c0) / *beta); tau = ei_conj((beta - c0) / beta);
} }
} }

2
Eigen/src/Householder/HouseholderSequence.h
View File

@ -162,7 +162,7 @@ private:
}; };
template<typename VectorsType, typename CoeffsType> template<typename VectorsType, typename CoeffsType>
HouseholderSequence<VectorsType,CoeffsType> makeHouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans=false) HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h, bool trans=false)
{ {
return HouseholderSequence<VectorsType,CoeffsType>(v, h, trans); return HouseholderSequence<VectorsType,CoeffsType>(v, h, trans);
} }

4
Eigen/src/QR/ColPivHouseholderQR.h
View File

@ -305,7 +305,7 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
} }
RealScalar beta; RealScalar beta;
m_qr.col(k).end(rows-k).makeHouseholderInPlace(&m_hCoeffs.coeffRef(k), &beta); m_qr.col(k).end(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
m_qr.coeffRef(k,k) = beta; m_qr.coeffRef(k,k) = beta;
m_qr.corner(BottomRight, rows-k, cols-k-1) m_qr.corner(BottomRight, rows-k, cols-k-1)
@ -347,7 +347,7 @@ struct ei_solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
typename Rhs::PlainMatrixType c(rhs()); typename Rhs::PlainMatrixType c(rhs());
// Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
c.applyOnTheLeft(makeHouseholderSequence( c.applyOnTheLeft(householderSequence(
dec().matrixQR().corner(TopLeft,rows,dec().rank()), dec().matrixQR().corner(TopLeft,rows,dec().rank()),
dec().hCoeffs().start(dec().rank())).transpose() dec().hCoeffs().start(dec().rank())).transpose()
); );

2
Eigen/src/QR/FullPivHouseholderQR.h
View File

@ -315,7 +315,7 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons
} }
RealScalar beta; RealScalar beta;
m_qr.col(k).end(rows-k).makeHouseholderInPlace(&m_hCoeffs.coeffRef(k), &beta); m_qr.col(k).end(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
m_qr.coeffRef(k,k) = beta; m_qr.coeffRef(k,k) = beta;
m_qr.corner(BottomRight, rows-k, cols-k-1) m_qr.corner(BottomRight, rows-k, cols-k-1)

4
Eigen/src/QR/HouseholderQR.h
View File

@ -198,7 +198,7 @@ HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType&
int remainingCols = cols - k - 1; int remainingCols = cols - k - 1;
RealScalar beta; RealScalar beta;
m_qr.col(k).end(remainingRows).makeHouseholderInPlace(&m_hCoeffs.coeffRef(k), &beta); m_qr.col(k).end(remainingRows).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
m_qr.coeffRef(k,k) = beta; m_qr.coeffRef(k,k) = beta;
// apply H to remaining part of m_qr from the left // apply H to remaining part of m_qr from the left
@ -225,7 +225,7 @@ struct ei_solve_retval<HouseholderQR<_MatrixType>, Rhs>
typename Rhs::PlainMatrixType c(rhs()); typename Rhs::PlainMatrixType c(rhs());
// Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
c.applyOnTheLeft(makeHouseholderSequence( c.applyOnTheLeft(householderSequence(
dec().matrixQR().corner(TopLeft,rows,rank), dec().matrixQR().corner(TopLeft,rows,rank),
dec().hCoeffs().start(rank)).transpose() dec().hCoeffs().start(rank)).transpose()
); );

6
test/householder.cpp
View File

@ -50,7 +50,7 @@ template<typename MatrixType> void householder(const MatrixType& m)
VectorType v1 = VectorType::Random(rows), v2; VectorType v1 = VectorType::Random(rows), v2;
v2 = v1; v2 = v1;
v1.makeHouseholder(&essential, &beta, &alpha); v1.makeHouseholder(essential, beta, alpha);
v1.applyHouseholderOnTheLeft(essential,beta,tmp); v1.applyHouseholderOnTheLeft(essential,beta,tmp);
VERIFY_IS_APPROX(v1.norm(), v2.norm()); VERIFY_IS_APPROX(v1.norm(), v2.norm());
VERIFY_IS_MUCH_SMALLER_THAN(v1.end(rows-1).norm(), v1.norm()); VERIFY_IS_MUCH_SMALLER_THAN(v1.end(rows-1).norm(), v1.norm());
@ -66,7 +66,7 @@ template<typename MatrixType> void householder(const MatrixType& m)
if(even) v1.end(rows-1).setZero(); if(even) v1.end(rows-1).setZero();
m1.colwise() = v1; m1.colwise() = v1;
m2 = m1; m2 = m1;
m1.col(0).makeHouseholder(&essential, &beta, &alpha); m1.col(0).makeHouseholder(essential, beta, alpha);
m1.applyHouseholderOnTheLeft(essential,beta,tmp); m1.applyHouseholderOnTheLeft(essential,beta,tmp);
VERIFY_IS_APPROX(m1.norm(), m2.norm()); VERIFY_IS_APPROX(m1.norm(), m2.norm());
VERIFY_IS_MUCH_SMALLER_THAN(m1.block(1,0,rows-1,cols).norm(), m1.norm()); VERIFY_IS_MUCH_SMALLER_THAN(m1.block(1,0,rows-1,cols).norm(), m1.norm());
@ -78,7 +78,7 @@ template<typename MatrixType> void householder(const MatrixType& m)
SquareMatrixType m3(rows,rows), m4(rows,rows); SquareMatrixType m3(rows,rows), m4(rows,rows);
m3.rowwise() = v1.transpose(); m3.rowwise() = v1.transpose();
m4 = m3; m4 = m3;
m3.row(0).makeHouseholder(&essential, &beta, &alpha); m3.row(0).makeHouseholder(essential, beta, alpha);
m3.applyHouseholderOnTheRight(essential,beta,tmp); m3.applyHouseholderOnTheRight(essential,beta,tmp);
VERIFY_IS_APPROX(m3.norm(), m4.norm()); VERIFY_IS_APPROX(m3.norm(), m4.norm());
VERIFY_IS_MUCH_SMALLER_THAN(m3.block(0,1,rows,rows-1).norm(), m3.norm()); VERIFY_IS_MUCH_SMALLER_THAN(m3.block(0,1,rows,rows-1).norm(), m3.norm());