devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix bug #314:

Browse Source 
- remove most of the metaprogramming kung fu in MathFunctions.h (only keep functions that differs from the std)
- remove the overloads for array expression that were in the std namespace
This commit is contained in:
Gael Guennebaud
2012-11-06 15:25:50 +01:00
parent 959ef37006
commit a76fbbf397
88 changed files with 496 additions and 468 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
Eigen/src/Eigenvalues/ComplexSchur.h
View File

@@ -258,10 +258,11 @@ inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i)
template<typename MatrixType>
typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::computeShift(Index iu, Index iter)
{
using std::abs;
if (iter == 10 || iter == 20)
{
// exceptional shift, taken from http://www.netlib.org/eispack/comqr.f
return internal::abs(internal::real(m_matT.coeff(iu,iu-1))) + internal::abs(internal::real(m_matT.coeff(iu-1,iu-2)));
return abs(internal::real(m_matT.coeff(iu,iu-1))) + abs(internal::real(m_matT.coeff(iu-1,iu-2)));
}
// compute the shift as one of the eigenvalues of t, the 2x2

20
Eigen/src/Eigenvalues/EigenSolver.h
View File

@@ -364,6 +364,8 @@ template<typename MatrixType>
EigenSolver<MatrixType>&
EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
{
using std::sqrt;
using std::abs;
assert(matrix.cols() == matrix.rows());
// Reduce to real Schur form.
@@ -388,7 +390,7 @@ EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvect
else
{
Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
Scalar z = internal::sqrt(internal::abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
Scalar z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
i += 2;
@@ -410,8 +412,9 @@ EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvect
template<typename Scalar>
std::complex<Scalar> cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi)
{
using std::abs;
Scalar r,d;
if (internal::abs(yr) > internal::abs(yi))
if (abs(yr) > abs(yi))
{
r = yi/yr;
d = yr + r*yi;
@@ -429,6 +432,7 @@ std::complex<Scalar> cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi)
template<typename MatrixType>
void EigenSolver<MatrixType>::doComputeEigenvectors()
{
using std::abs;
const Index size = m_eivec.cols();
const Scalar eps = NumTraits<Scalar>::epsilon();
@@ -484,14 +488,14 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
Scalar t = (x * lastr - lastw * r) / denom;
m_matT.coeffRef(i,n) = t;
if (internal::abs(x) > internal::abs(lastw))
if (abs(x) > abs(lastw))
m_matT.coeffRef(i+1,n) = (-r - w * t) / x;
else
m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw;
}
// Overflow control
Scalar t = internal::abs(m_matT.coeff(i,n));
Scalar t = abs(m_matT.coeff(i,n));
if ((eps * t) * t > Scalar(1))
m_matT.col(n).tail(size-i) /= t;
}
@@ -503,7 +507,7 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
Index l = n-1;
// Last vector component imaginary so matrix is triangular
if (internal::abs(m_matT.coeff(n,n-1)) > internal::abs(m_matT.coeff(n-1,n)))
if (abs(m_matT.coeff(n,n-1)) > abs(m_matT.coeff(n-1,n)))
{
m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1);
m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1);
@@ -545,12 +549,12 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
Scalar vr = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag() - q * q;
Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
if ((vr == 0.0) && (vi == 0.0))
vr = eps * norm * (internal::abs(w) + internal::abs(q) + internal::abs(x) + internal::abs(y) + internal::abs(lastw));
vr = eps * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(lastw));
std::complex<Scalar> cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi);
m_matT.coeffRef(i,n-1) = internal::real(cc);
m_matT.coeffRef(i,n) = internal::imag(cc);
if (internal::abs(x) > (internal::abs(lastw) + internal::abs(q)))
if (abs(x) > (abs(lastw) + abs(q)))
{
m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x;
m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x;
@@ -565,7 +569,7 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
// Overflow control
using std::max;
Scalar t = (max)(internal::abs(m_matT.coeff(i,n-1)),internal::abs(m_matT.coeff(i,n)));
Scalar t = (max)(abs(m_matT.coeff(i,n-1)),abs(m_matT.coeff(i,n)));
if ((eps * t) * t > Scalar(1))
m_matT.block(i, n-1, size-i, 2) /= t;

4
Eigen/src/Eigenvalues/GeneralizedEigenSolver.h
View File

@@ -290,6 +290,8 @@ template<typename MatrixType>
GeneralizedEigenSolver<MatrixType>&
GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors)
{
using std::sqrt;
using std::abs;
eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows());
// Reduce to generalized real Schur form:
@@ -317,7 +319,7 @@ GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixTyp
else
{
Scalar p = Scalar(0.5) * (m_matS.coeff(i, i) - m_matS.coeff(i+1, i+1));
Scalar z = internal::sqrt(internal::abs(p * p + m_matS.coeff(i+1, i) * m_matS.coeff(i, i+1)));
Scalar z = sqrt(abs(p * p + m_matS.coeff(i+1, i) * m_matS.coeff(i, i+1)));
m_alphas.coeffRef(i) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, z);
m_alphas.coeffRef(i+1) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, -z);

3
Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h
View File

@@ -121,10 +121,11 @@ template<typename Derived>
inline typename MatrixBase<Derived>::RealScalar
MatrixBase<Derived>::operatorNorm() const
{
using std::sqrt;
typename Derived::PlainObject m_eval(derived());
// FIXME if it is really guaranteed that the eigenvalues are already sorted,
// then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough.
return internal::sqrt((m_eval*m_eval.adjoint())
return sqrt((m_eval*m_eval.adjoint())
.eval()
.template selfadjointView<Lower>()
.eigenvalues()

20
Eigen/src/Eigenvalues/RealQZ.h
View File

@@ -278,13 +278,14 @@ namespace Eigen {
template<typename MatrixType>
inline typename MatrixType::Index RealQZ<MatrixType>::findSmallSubdiagEntry(Index iu)
{
using std::abs;
Index res = iu;
while (res > 0)
{
Scalar s = internal::abs(m_S.coeff(res-1,res-1)) + internal::abs(m_S.coeff(res,res));
Scalar s = abs(m_S.coeff(res-1,res-1)) + abs(m_S.coeff(res,res));
if (s == Scalar(0.0))
s = m_normOfS;
if (internal::abs(m_S.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
if (abs(m_S.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
break;
res--;
}
@@ -295,9 +296,10 @@ namespace Eigen {
template<typename MatrixType>
inline typename MatrixType::Index RealQZ<MatrixType>::findSmallDiagEntry(Index f, Index l)
{
using std::abs;
Index res = l;
while (res >= f) {
if (internal::abs(m_T.coeff(res,res)) <= NumTraits<Scalar>::epsilon() * m_normOfT)
if (abs(m_T.coeff(res,res)) <= NumTraits<Scalar>::epsilon() * m_normOfT)
break;
res--;
}
@@ -308,8 +310,10 @@ namespace Eigen {
template<typename MatrixType>
inline void RealQZ<MatrixType>::splitOffTwoRows(Index i)
{
using std::abs;
using std::sqrt;
const Index dim=m_S.cols();
if (internal::abs(m_S.coeff(i+1,i)==Scalar(0)))
if (abs(m_S.coeff(i+1,i)==Scalar(0)))
return;
Index z = findSmallDiagEntry(i,i+1);
if (z==i-1)
@@ -320,7 +324,7 @@ namespace Eigen {
Scalar p = Scalar(0.5)*(STi(0,0)-STi(1,1));
Scalar q = p*p + STi(1,0)*STi(0,1);
if (q>=0) {
Scalar z = internal::sqrt(q);
Scalar z = sqrt(q);
// one QR-like iteration for ABi - lambda I
// is enough - when we know exact eigenvalue in advance,
// convergence is immediate
@@ -393,7 +397,9 @@ namespace Eigen {
/** \internal QR-like iterative step for block f..l */
template<typename MatrixType>
inline void RealQZ<MatrixType>::step(Index f, Index l, Index iter) {
inline void RealQZ<MatrixType>::step(Index f, Index l, Index iter)
{
using std::abs;
const Index dim = m_S.cols();
// x, y, z
@@ -411,7 +417,7 @@ namespace Eigen {
a98=m_S.coeff(l-0,l-1),
b77i=Scalar(1.0)/m_T.coeff(l-2,l-2),
b88i=Scalar(1.0)/m_T.coeff(l-1,l-1);
Scalar ss = internal::abs(a87*b77i) + internal::abs(a98*b88i),
Scalar ss = abs(a87*b77i) + abs(a98*b88i),
lpl = Scalar(1.5)*ss,
ll = ss*ss;
x = ll + a11*a11*b11i*b11i - lpl*a11*b11i + a12*a21*b11i*b22i

22
Eigen/src/Eigenvalues/RealSchur.h
View File

@@ -345,13 +345,14 @@ inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT()
template<typename MatrixType>
inline typename MatrixType::Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu, Scalar norm)
{
using std::abs;
Index res = iu;
while (res > 0)
{
Scalar s = internal::abs(m_matT.coeff(res-1,res-1)) + internal::abs(m_matT.coeff(res,res));
Scalar s = abs(m_matT.coeff(res-1,res-1)) + abs(m_matT.coeff(res,res));
if (s == 0.0)
s = norm;
if (internal::abs(m_matT.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
if (abs(m_matT.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
break;
res--;
}
@@ -362,6 +363,8 @@ inline typename MatrixType::Index RealSchur<MatrixType>::findSmallSubdiagEntry(I
template<typename MatrixType>
inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, Scalar exshift)
{
using std::sqrt;
using std::abs;
const Index size = m_matT.cols();
// The eigenvalues of the 2x2 matrix [a b; c d] are
@@ -373,7 +376,7 @@ inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, Scal
if (q >= Scalar(0)) // Two real eigenvalues
{
Scalar z = internal::sqrt(internal::abs(q));
Scalar z = sqrt(abs(q));
JacobiRotation<Scalar> rot;
if (p >= Scalar(0))
rot.makeGivens(p + z, m_matT.coeff(iu, iu-1));
@@ -395,6 +398,8 @@ inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, Scal
template<typename MatrixType>
inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo)
{
using std::sqrt;
using std::abs;
shiftInfo.coeffRef(0) = m_matT.coeff(iu,iu);
shiftInfo.coeffRef(1) = m_matT.coeff(iu-1,iu-1);
shiftInfo.coeffRef(2) = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
@@ -405,7 +410,7 @@ inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& ex
exshift += shiftInfo.coeff(0);
for (Index i = 0; i <= iu; ++i)
m_matT.coeffRef(i,i) -= shiftInfo.coeff(0);
Scalar s = internal::abs(m_matT.coeff(iu,iu-1)) + internal::abs(m_matT.coeff(iu-1,iu-2));
Scalar s = abs(m_matT.coeff(iu,iu-1)) + abs(m_matT.coeff(iu-1,iu-2));
shiftInfo.coeffRef(0) = Scalar(0.75) * s;
shiftInfo.coeffRef(1) = Scalar(0.75) * s;
shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s;
@@ -418,7 +423,7 @@ inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& ex
s = s * s + shiftInfo.coeff(2);
if (s > Scalar(0))
{
s = internal::sqrt(s);
s = sqrt(s);
if (shiftInfo.coeff(1) < shiftInfo.coeff(0))
s = -s;
s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
@@ -435,6 +440,7 @@ inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& ex
template<typename MatrixType>
inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector)
{
using std::abs;
Vector3s& v = firstHouseholderVector; // alias to save typing
for (im = iu-2; im >= il; --im)
@@ -448,9 +454,9 @@ inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const V
if (im == il) {
break;
}
const Scalar lhs = m_matT.coeff(im,im-1) * (internal::abs(v.coeff(1)) + internal::abs(v.coeff(2)));
const Scalar rhs = v.coeff(0) * (internal::abs(m_matT.coeff(im-1,im-1)) + internal::abs(Tmm) + internal::abs(m_matT.coeff(im+1,im+1)));
if (internal::abs(lhs) < NumTraits<Scalar>::epsilon() * rhs)
const Scalar lhs = m_matT.coeff(im,im-1) * (abs(v.coeff(1)) + abs(v.coeff(2)));
const Scalar rhs = v.coeff(0) * (abs(m_matT.coeff(im-1,im-1)) + abs(Tmm) + abs(m_matT.coeff(im+1,im+1)));
if (abs(lhs) < NumTraits<Scalar>::epsilon() * rhs)
{
break;
}

5
Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
View File

@@ -384,6 +384,7 @@ template<typename MatrixType>
SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
::compute(const MatrixType& matrix, int options)
{
using std::abs;
eigen_assert(matrix.cols() == matrix.rows());
eigen_assert((options&~(EigVecMask|GenEigMask))==0
&& (options&EigVecMask)!=EigVecMask
@@ -421,7 +422,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
while (end>0)
{
for (Index i = start; i<end; ++i)
if (internal::isMuchSmallerThan(internal::abs(m_subdiag[i]),(internal::abs(diag[i])+internal::abs(diag[i+1]))))
if (internal::isMuchSmallerThan(abs(m_subdiag[i]),(abs(diag[i])+abs(diag[i+1]))))
m_subdiag[i] = 0;
// find the largest unreduced block
@@ -675,6 +676,7 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2
static inline void run(SolverType& solver, const MatrixType& mat, int options)
{
using std::sqrt;
eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows());
eigen_assert((options&~(EigVecMask|GenEigMask))==0
&& (options&EigVecMask)!=EigVecMask
@@ -736,6 +738,7 @@ namespace internal {
template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n)
{
using std::abs;
RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
RealScalar e = subdiag[end-1];
// Note that thanks to scaling, e^2 or td^2 cannot overflow, however they can still

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

@@ -345,6 +345,7 @@ namespace internal {
template<typename MatrixType, typename CoeffVectorType>
void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs)
{
using internal::conj;
typedef typename MatrixType::Index Index;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
@@ -467,6 +468,7 @@ struct tridiagonalization_inplace_selector<MatrixType,3,false>
template<typename DiagonalType, typename SubDiagonalType>
static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
{
using std::sqrt;
diag[0] = mat(0,0);
RealScalar v1norm2 = abs2(mat(2,0));
if(v1norm2 == RealScalar(0))