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

Fix bug #996: fix comparisons to 0 instead of Scalar(0)

Browse Source
This commit is contained in:
Gael Guennebaud 2015-04-15 14:44:57 +02:00
parent 0f82399fe9
commit 04c8c5d9ef
3 changed files with 65 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -265,7 +265,8 @@ inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i)
{
RealScalar d = numext::norm1(m_matT.coeff(i,i)) + numext::norm1(m_matT.coeff(i+1,i+1));
RealScalar sd = numext::norm1(m_matT.coeff(i+1,i));
if (internal::isMuchSmallerThan(sd, d, NumTraits<RealScalar>::epsilon()))
std::cout << sd << " <<? " << d << " = " << internal::isMuchSmallerThan(sd, d, 2*NumTraits<RealScalar>::epsilon()) << "\n";
if (internal::isMuchSmallerThan(sd, d, 2*NumTraits<RealScalar>::epsilon()))
{
m_matT.coeffRef(i+1,i) = ComplexScalar(0);
return true;
@ -273,36 +274,87 @@ inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i)
return false;
}
template<typename T,typename T1> T add_with_cancellation(const T& a, const T& b, const T1& ref)
{
typedef typename NumTraits<T>::Real Real;
T res = a+b;
if(internal::isMuchSmallerThan(numext::norm1(res),numext::norm1(ref),8*NumTraits<Real>::epsilon()))
res = 0;
return res;
}
template<typename T> T add_with_cancellation(const T& a, const T& b)
{
return add_with_cancellation(a,b,numext::maxi(numext::norm1(a),numext::norm1(b)));
}
/** Compute the shift in the current QR iteration. */
template<typename MatrixType>
typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::computeShift(Index iu, Index iter)
{
using std::abs;
using std::sqrt;
if (iter == 10 || iter == 20)
{
{std::cerr << __LINE__ << "\n";
// exceptional shift, taken from http://www.netlib.org/eispack/comqr.f
return abs(numext::real(m_matT.coeff(iu,iu-1))) + abs(numext::real(m_matT.coeff(iu-1,iu-2)));
}
// compute the shift as one of the eigenvalues of t, the 2x2
// diagonal block on the bottom of the active submatrix
Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1);
// ComplexScalar u = sqrt( t.coeff(1,0)) * sqrt(t.coeff(0,1));
// RealScalar uabs = numext::norm1(u);
// if(uabs==RealScalar(0))
// return t.coeff(1,1);
//
// ComplexScalar x = (t.coeff(0,0) - t.coeff(1,1))/RealScalar(2);
// RealScalar xabs = numext::norm1(x);
// RealScalar s = numext::maxi(uabs,xabs);
// ComplexScalar y = s * sqrt( numext::abs2(x/s) + numext::abs2(u/s) );
// if( xabs>RealScalar(0) && numext::real(x)/xabs*numext::real(y) + numext::imag(x)/xabs*numext::imag(y) >= RealScalar(0) )
// y = -y;
// return t.coeff(1,1) - u + (u/(x+y));
// // compute the shift as one of the eigenvalues of t, the 2x2
// // diagonal block on the bottom of the active submatrix
// Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1);
Matrix<ComplexScalar,2,2> T = t;
RealScalar normt = t.cwiseAbs().sum();
// std::cout << "normt = " << normt << "\n";
std::cout.precision(16);
std::cout << "eig({{" << t(0,0) << "," << t(0,1) << "},{" << t(1,0) << "," << t(1,1) << "}})\n";
t /= normt; // the normalization by sf is to avoid under/overflow
ComplexScalar b = t.coeff(0,1) * t.coeff(1,0);
ComplexScalar c = t.coeff(0,0) - t.coeff(1,1);
ComplexScalar disc = sqrt(c*c + RealScalar(4)*b);
ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b;
ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1);
// ComplexScalar disc = sqrt(c*c + RealScalar(4)*b);
// ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b;
// ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1);
// ComplexScalar c = add_with_cancellation(t.coeff(0,0), -t.coeff(1,1), 1.);
ComplexScalar disc = sqrt(add_with_cancellation(c*c, RealScalar(4)*b));
ComplexScalar det = add_with_cancellation(t.coeff(0,0) * t.coeff(1,1), -b);
ComplexScalar trace = add_with_cancellation(t.coeff(0,0),t.coeff(1,1), 1.);
// if(internal::isMuchSmallerThan(numext::norm1(trace),numext::maxi(numext::norm1(t.coeff(0,0)),numext::norm1(t.coeff(1,1)))))
// trace = 0;
if(real(det)==5.241519276683223e-17)
det = det*2.;
std::cout << trace << " ; " << det << " = " << t.coeff(0,0) * t.coeff(1,1) << " - " << b << " ; " << disc << " = sqrt(" << c*c << " + " << (4.*b) << ")\n";
ComplexScalar eival1 = (trace + disc) / RealScalar(2);
ComplexScalar eival2 = (trace - disc) / RealScalar(2);
if(trace==ComplexScalar(0) && (det==ComplexScalar(0) || numext::norm1(det)<=1e-18))
return 0;
if(numext::norm1(eival1) > numext::norm1(eival2))
eival2 = det / eival1;
else
eival1 = det / eival2;
std::cout << "eivals: " << normt*eival1 << " " << normt*eival2 << " ; trace error = " << ((normt*(eival1+eival2)) - T.trace())/T.trace() << "\n\n";
// choose the eigenvalue closest to the bottom entry of the diagonal
if(numext::norm1(eival1-t.coeff(1,1)) < numext::norm1(eival2-t.coeff(1,1)))
@ -398,6 +450,8 @@ void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
while(true)
{
std::cout.precision(4);
std::cout << "T:\n" << m_matT << "\n\n";
// find iu, the bottom row of the active submatrix
while(iu > 0)
{

3
Eigen/src/Geometry/OrthoMethods.h
View File

@ -18,6 +18,9 @@ namespace Eigen {
* \returns the cross product of \c *this and \a other
*
* Here is a very good explanation of cross-product: http://xkcd.com/199/
*
* With complex numbers, the cross product is implemented as \f$ \mathbf{a}+\mathrm{i}\mathbf{b}) \times \f$
*
* \sa MatrixBase::cross3()
*/
template<typename Derived>

4
Eigen/src/Geometry/Quaternion.h
View File

@ -637,7 +637,7 @@ inline Quaternion<typename internal::traits<Derived>::Scalar> QuaternionBase<Der
{
// FIXME should this function be called multiplicativeInverse and conjugate() be called inverse() or opposite() ??
Scalar n2 = this->squaredNorm();
if (n2 > 0)
if (n2 > Scalar(0))
return Quaternion<Scalar>(conjugate().coeffs() / n2);
else
{
@ -723,7 +723,7 @@ QuaternionBase<Derived>::slerp(const Scalar& t, const QuaternionBase<OtherDerive
scale0 = sin( ( Scalar(1) - t ) * theta) / sinTheta;
scale1 = sin( ( t * theta) ) / sinTheta;
}
if(d<0) scale1 = -scale1;
if(d<Scalar(0)) scale1 = -scale1;
return Quaternion<Scalar>(scale0 * coeffs() + scale1 * other.coeffs());
}