diff --git a/Eigen/src/Jacobi/Jacobi.h b/Eigen/src/Jacobi/Jacobi.h index c38d4583f..66ba06b19 100644 --- a/Eigen/src/Jacobi/Jacobi.h +++ b/Eigen/src/Jacobi/Jacobi.h @@ -195,7 +195,7 @@ void PlanarRotation::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$ diff --git a/Eigen/src/QR/ComplexEigenSolver.h b/Eigen/src/QR/ComplexEigenSolver.h index af47c2195..6caa6bc2d 100644 --- a/Eigen/src/QR/ComplexEigenSolver.h +++ b/Eigen/src/QR/ComplexEigenSolver.h @@ -107,7 +107,7 @@ void ComplexEigenSolver::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; } diff --git a/Eigen/src/QR/ComplexSchur.h b/Eigen/src/QR/ComplexSchur.h index 153826811..b60be5430 100644 --- a/Eigen/src/QR/ComplexSchur.h +++ b/Eigen/src/QR/ComplexSchur.h @@ -80,6 +80,43 @@ template class ComplexSchur bool m_isInitialized; }; +/** Computes the principal value of the square root of the complex \a z. */ +template +std::complex ei_sqrt(const std::complex &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(tre,tim)); + +} + template void ComplexSchur::compute(const MatrixType& matrix) { @@ -146,7 +183,7 @@ void ComplexSchur::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::compute(const MatrixType& matrix) } // FIXME : is it necessary ? + /* for(int i=0 ; i -std::complex csqrt(const std::complex &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(tre,tim)); - -} - - #endif // EIGEN_COMPLEX_SCHUR_H