fix #51 (bad use of std::complex::real)

3 files changed, 46 insertions, 47 deletions

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<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$
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<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;

 

     }

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<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));
 
+  m_isInitialized = true;
 }
 
-
 #endif // EIGEN_COMPLEX_SCHUR_H