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 #51 (bad use of std::complex::real)

Browse Source
This commit is contained in:
Gael Guennebaud 2009-09-02 15:18:11 +02:00
parent b83654b5d0
commit 7586f7f706
3 changed files with 46 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Eigen/src/Jacobi/Jacobi.h
View File

@ -195,7 +195,7 @@ void PlanarRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar
/****************************************************************************************
* Implementation of MatrixBase methods
/***************************************************************************************/
****************************************************************************************/
/** Applies the clock wise 2D rotation \a j to the set of 2D vectors of cordinates \a x and \a y:
* \f$ \left ( \begin{array}{cc} x \\ y \end{array} \right ) = J \left ( \begin{array}{cc} x \\ y \end{array} \right ) \f$

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

@ -107,7 +107,7 @@ void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
m_eivec.coeffRef(i,k) -= (schur.matrixT().row(i).segment(i+1,k-i-1) * m_eivec.col(k).segment(i+1,k-i-1)).value();
z = schur.matrixT().coeff(i,i) - d2;
if(z==Scalar(0))
z.real() = eps * norm;
ei_real_ref(z) = eps * norm;
m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) / z;
}

89
Eigen/src/QR/ComplexSchur.h
View File

@ -80,6 +80,43 @@ template<typename _MatrixType> class ComplexSchur
bool m_isInitialized;
};
/** Computes the principal value of the square root of the complex \a z. */
template<typename RealScalar>
std::complex<RealScalar> ei_sqrt(const std::complex<RealScalar> &z)
{
RealScalar t, tre, tim;
t = ei_abs(z);
if (ei_abs(ei_real(z)) <= ei_abs(ei_real(z)))
{
// No cancellation in these formulas
tre = ei_sqrt(0.5*(t + ei_real(z)));
tim = ei_sqrt(0.5*(t - ei_real(z)));
}
else
{
// Stable computation of the above formulas
if (z.real() > 0)
{
tre = t + z.real();
tim = ei_abs(ei_imag(z))*ei_sqrt(0.5/tre);
tre = ei_sqrt(0.5*tre);
}
else
{
tim = t - z.real();
tre = ei_abs(ei_imag(z))*ei_sqrt(0.5/tim);
tim = ei_sqrt(0.5*tim);
}
}
if(z.imag() < 0)
tim = -tim;
return (std::complex<RealScalar>(tre,tim));
}
template<typename MatrixType>
void ComplexSchur<MatrixType>::compute(const MatrixType& matrix)
{
@ -146,7 +183,7 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix)
c = t.determinant();
b = t.diagonal().sum();
disc = csqrt(b*b - RealScalar(4)*c);
disc = ei_sqrt(b*b - RealScalar(4)*c);
r1 = (b+disc)/RealScalar(2);
r2 = (b-disc)/RealScalar(2);
@ -183,56 +220,18 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix)
}
// FIXME : is it necessary ?
/*
for(int i=0 ; i<n ; i++)
for(int j=0 ; j<n ; j++)
{
if(ei_abs(m_matT.coeff(i,j).real()) < eps)
m_matT.coeffRef(i,j).real() = 0;
if(ei_abs(m_matT.coeff(i,j).imag()) < eps)
m_matT.coeffRef(i,j).imag() = 0;
if(ei_abs(ei_real(m_matT.coeff(i,j))) < eps)
ei_real_ref(m_matT.coeffRef(i,j)) = 0;
if(ei_imag(ei_abs(m_matT.coeff(i,j))) < eps)
ei_imag_ref(m_matT.coeffRef(i,j)) = 0;
}
*/
m_isInitialized = true;
}
/**
* Computes the principal value of the square root of the complex \a z.
*/
template<typename RealScalar>
std::complex<RealScalar> csqrt(const std::complex<RealScalar> &z)
{
RealScalar t, tre, tim;
t = ei_abs(z);
if (ei_abs(z.real()) <= ei_abs(z.imag()))
{
// No cancellation in these formulas
tre = ei_sqrt(0.5*(t+z.real()));
tim = ei_sqrt(0.5*(t-z.real()));
}
else
{
// Stable computation of the above formulas
if (z.real() > 0)
{
tre = t + z.real();
tim = ei_abs(z.imag())*ei_sqrt(0.5/tre);
tre = ei_sqrt(0.5*tre);
}
else
{
tim = t - z.real();
tre = ei_abs(z.imag())*ei_sqrt(0.5/tim);
tim = ei_sqrt(0.5*tim);
}
}
if(z.imag() < 0)
tim = -tim;
return (std::complex<RealScalar>(tre,tim));
}
#endif // EIGEN_COMPLEX_SCHUR_H