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

more Scalar conversions fixes

Browse Source
This commit is contained in:
Gael Guennebaud 2009-01-07 10:22:46 +00:00
parent ac9d805c2f
commit 94efa16187
3 changed files with 16 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Eigen/src/Core/MathFunctions.h
View File

@ -94,6 +94,7 @@ inline float ei_log(float x) { return std::log(x); }
inline float ei_sin(float x) { return std::sin(x); } inline float ei_sin(float x) { return std::sin(x); }
inline float ei_cos(float x) { return std::cos(x); } inline float ei_cos(float x) { return std::cos(x); }
inline float ei_pow(float x, float y) { return std::pow(x, y); } inline float ei_pow(float x, float y) { return std::pow(x, y); }
inline float ei_hypot(float x, float y) { return hypotf(x,y); }
template<> inline float ei_random(float a, float b) template<> inline float ei_random(float a, float b)
{ {
@ -139,6 +140,7 @@ inline double ei_log(double x) { return std::log(x); }
inline double ei_sin(double x) { return std::sin(x); } inline double ei_sin(double x) { return std::sin(x); }
inline double ei_cos(double x) { return std::cos(x); } inline double ei_cos(double x) { return std::cos(x); }
inline double ei_pow(double x, double y) { return std::pow(x, y); } inline double ei_pow(double x, double y) { return std::pow(x, y); }
inline double ei_hypot(double x, double y) { return hypot(x,y); }
template<> inline double ei_random(double a, double b) template<> inline double ei_random(double a, double b)
{ {

12
Eigen/src/QR/EigenSolver.h
View File

@ -282,7 +282,7 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
int n = nn-1; int n = nn-1;
int low = 0; int low = 0;
int high = nn-1; int high = nn-1;
Scalar eps = Scalar(pow(2.0,-52.0)); Scalar eps = ei_pow(Scalar(2),ei_is_same_type<Scalar,float>::ret ? Scalar(-23) : Scalar(-52));
Scalar exshift = 0.0; Scalar exshift = 0.0;
Scalar p=0,q=0,r=0,s=0,z=0,t,w,x,y; Scalar p=0,q=0,r=0,s=0,z=0,t,w,x,y;
@ -328,7 +328,7 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
else if (l == n-1) // Two roots found else if (l == n-1) // Two roots found
{ {
w = matH.coeff(n,n-1) * matH.coeff(n-1,n); w = matH.coeff(n,n-1) * matH.coeff(n-1,n);
p = Scalar((matH.coeff(n-1,n-1) - matH.coeff(n,n)) / 2.0); p = (matH.coeff(n-1,n-1) - matH.coeff(n,n)) * Scalar(0.5);
q = p * p + w; q = p * p + w;
z = ei_sqrt(ei_abs(q)); z = ei_sqrt(ei_abs(q));
matH.coeffRef(n,n) = matH.coeff(n,n) + exshift; matH.coeffRef(n,n) = matH.coeff(n,n) + exshift;
@ -405,8 +405,8 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
for (int i = low; i <= n; ++i) for (int i = low; i <= n; ++i)
matH.coeffRef(i,i) -= x; matH.coeffRef(i,i) -= x;
s = ei_abs(matH.coeff(n,n-1)) + ei_abs(matH.coeff(n-1,n-2)); s = ei_abs(matH.coeff(n,n-1)) + ei_abs(matH.coeff(n-1,n-2));
x = y = Scalar(0.75 * s); x = y = Scalar(0.75) * s;
w = Scalar(-0.4375 * s * s); w = Scalar(-0.4375) * s * s;
} }
// MATLAB's new ad hoc shift // MATLAB's new ad hoc shift
@ -469,7 +469,7 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
if (k != m) { if (k != m) {
p = matH.coeff(k,k-1); p = matH.coeff(k,k-1);
q = matH.coeff(k+1,k-1); q = matH.coeff(k+1,k-1);
r = Scalar(notlast ? matH.coeff(k+2,k-1) : 0.0); r = notlast ? matH.coeff(k+2,k-1) : Scalar(0);
x = ei_abs(p) + ei_abs(q) + ei_abs(r); x = ei_abs(p) + ei_abs(q) + ei_abs(r);
if (x != 0.0) if (x != 0.0)
{ {
@ -647,7 +647,7 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
x = matH.coeff(i,i+1); x = matH.coeff(i,i+1);
y = matH.coeff(i+1,i); y = matH.coeff(i+1,i);
vr = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag() - q * q; vr = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag() - q * q;
vi = Scalar((m_eivalues.coeff(i).real() - p) * 2.0 * q); vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
if ((vr == 0.0) && (vi == 0.0)) if ((vr == 0.0) && (vi == 0.0))
vr = eps * norm * (ei_abs(w) + ei_abs(q) + ei_abs(x) + ei_abs(y) + ei_abs(z)); vr = eps * norm * (ei_abs(w) + ei_abs(q) + ei_abs(x) + ei_abs(y) + ei_abs(z));

16
Eigen/src/SVD/SVD.h
View File

@ -242,7 +242,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
// Main iteration loop for the singular values. // Main iteration loop for the singular values.
int pp = p-1; int pp = p-1;
int iter = 0; int iter = 0;
Scalar eps(Scalar(pow(2.0,-52.0))); Scalar eps = ei_pow(Scalar(2),ei_is_same_type<Scalar,float>::ret ? Scalar(-23) : Scalar(-52));
while (p > 0) while (p > 0)
{ {
int k=0; int k=0;
@ -281,7 +281,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
{ {
if (ks == k) if (ks == k)
break; break;
Scalar t( Scalar((ks != p ? ei_abs(e[ks]) : 0.) + (ks != k+1 ? ei_abs(e[ks-1]) : 0.)) ); Scalar t = (ks != p ? ei_abs(e[ks]) : Scalar(0)) + (ks != k+1 ? ei_abs(e[ks-1]) : Scalar(0));
if (ei_abs(m_sigma[ks]) <= eps*t) if (ei_abs(m_sigma[ks]) <= eps*t)
{ {
m_sigma[ks] = 0.0; m_sigma[ks] = 0.0;
@ -315,7 +315,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
e[p-2] = 0.0; e[p-2] = 0.0;
for (j = p-2; j >= k; --j) for (j = p-2; j >= k; --j)
{ {
Scalar t(Scalar(hypot(m_sigma[j],f))); Scalar t(ei_hypot(m_sigma[j],f));
Scalar cs(m_sigma[j]/t); Scalar cs(m_sigma[j]/t);
Scalar sn(f/t); Scalar sn(f/t);
m_sigma[j] = t; m_sigma[j] = t;
@ -344,7 +344,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
e[k-1] = 0.0; e[k-1] = 0.0;
for (j = k; j < p; ++j) for (j = k; j < p; ++j)
{ {
Scalar t(Scalar(hypot(m_sigma[j],f))); Scalar t(ei_hypot(m_sigma[j],f));
Scalar cs( m_sigma[j]/t); Scalar cs( m_sigma[j]/t);
Scalar sn(f/t); Scalar sn(f/t);
m_sigma[j] = t; m_sigma[j] = t;
@ -375,7 +375,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
Scalar epm1 = e[p-2]/scale; Scalar epm1 = e[p-2]/scale;
Scalar sk = m_sigma[k]/scale; Scalar sk = m_sigma[k]/scale;
Scalar ek = e[k]/scale; Scalar ek = e[k]/scale;
Scalar b = Scalar(((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0); Scalar b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/Scalar(2);
Scalar c = (sp*epm1)*(sp*epm1); Scalar c = (sp*epm1)*(sp*epm1);
Scalar shift = 0.0; Scalar shift = 0.0;
if ((b != 0.0) || (c != 0.0)) if ((b != 0.0) || (c != 0.0))
@ -392,7 +392,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
for (j = k; j < p-1; ++j) for (j = k; j < p-1; ++j)
{ {
Scalar t = Scalar(hypot(f,g)); Scalar t = ei_hypot(f,g);
Scalar cs = f/t; Scalar cs = f/t;
Scalar sn = g/t; Scalar sn = g/t;
if (j != k) if (j != k)
@ -410,7 +410,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
m_matV(i,j) = t; m_matV(i,j) = t;
} }
} }
t = Scalar(hypot(f,g)); t = ei_hypot(f,g);
cs = f/t; cs = f/t;
sn = g/t; sn = g/t;
m_sigma[j] = t; m_sigma[j] = t;
@ -439,7 +439,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
// Make the singular values positive. // Make the singular values positive.
if (m_sigma[k] <= 0.0) if (m_sigma[k] <= 0.0)
{ {
m_sigma[k] = Scalar((m_sigma[k] < 0.0 ? -m_sigma[k] : 0.0)); m_sigma[k] = m_sigma[k] < Scalar(0) ? -m_sigma[k] : Scalar(0);
if (wantv) if (wantv)
m_matV.col(k).start(pp+1) = -m_matV.col(k).start(pp+1); m_matV.col(k).start(pp+1) = -m_matV.col(k).start(pp+1);
} }