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

Remove Eigen::internal::sqrt(), see bug #264.

Browse Source
This commit is contained in:
Jitse Niesen 2011-05-12 16:52:56 +01:00
parent 0aa7425f15
commit e22a523021
1 changed files with 1 additions and 41 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -227,46 +227,6 @@ template<typename _MatrixType> class ComplexSchur
friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>; friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>;
}; };
namespace internal {
/** Computes the principal value of the square root of the complex \a z. */
template<typename RealScalar>
std::complex<RealScalar> sqrt(const std::complex<RealScalar> &z)
{
RealScalar t, tre, tim;
t = abs(z);
if (abs(real(z)) <= abs(imag(z)))
{
// No cancellation in these formulas
tre = sqrt(RealScalar(0.5)*(t + real(z)));
tim = sqrt(RealScalar(0.5)*(t - real(z)));
}
else
{
// Stable computation of the above formulas
if (z.real() > RealScalar(0))
{
tre = t + z.real();
tim = abs(imag(z))*sqrt(RealScalar(0.5)/tre);
tre = sqrt(RealScalar(0.5)*tre);
}
else
{
tim = t - z.real();
tre = abs(imag(z))*sqrt(RealScalar(0.5)/tim);
tim = sqrt(RealScalar(0.5)*tim);
}
}
if(z.imag() < RealScalar(0))
tim = -tim;
return (std::complex<RealScalar>(tre,tim));
}
} // end namespace internal
/** If m_matT(i+1,i) is neglegible in floating point arithmetic /** If m_matT(i+1,i) is neglegible in floating point arithmetic
* compared to m_matT(i,i) and m_matT(j,j), then set it to zero and * compared to m_matT(i,i) and m_matT(j,j), then set it to zero and
* return true, else return false. */ * return true, else return false. */
@ -302,7 +262,7 @@ typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::compu
ComplexScalar b = t.coeff(0,1) * t.coeff(1,0); ComplexScalar b = t.coeff(0,1) * t.coeff(1,0);
ComplexScalar c = t.coeff(0,0) - t.coeff(1,1); ComplexScalar c = t.coeff(0,0) - t.coeff(1,1);
ComplexScalar disc = internal::sqrt(c*c + RealScalar(4)*b); ComplexScalar disc = sqrt(c*c + RealScalar(4)*b);
ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b; ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b;
ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1); ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1);
ComplexScalar eival1 = (trace + disc) / RealScalar(2); ComplexScalar eival1 = (trace + disc) / RealScalar(2);