mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-30 17:40:05 +08:00
Code Issues Packages Projects Releases Wiki Activity

updated EigenSolver to use .coeff / .coeffRef

Browse Source
This commit is contained in:
Gael Guennebaud 2008-05-31 16:31:10 +00:00
parent c9fb248c36
commit a2f71f9d7e
1 changed files with 161 additions and 161 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -162,7 +162,7 @@ void EigenSolver<MatrixType,IsSelfadjoint>::tridiagonalization(RealVectorType& e
if (scale == 0.0)
{
eivali[i] = eivalr[i-1];
eivali.coeffRef(i) = eivalr.coeff(i-1);
eivalr.start(i) = m_eivec.row(i-1).start(i);
m_eivec.corner(TopLeft, i, i) = m_eivec.corner(TopLeft, i, i).diagonal().asDiagonal();
}
@ -172,28 +172,28 @@ void EigenSolver<MatrixType,IsSelfadjoint>::tridiagonalization(RealVectorType& e
eivalr.start(i) /= scale;
h = eivalr.start(i).cwiseAbs2().sum();
Scalar f = eivalr[i-1];
Scalar f = eivalr.coeff(i-1);
Scalar g = ei_sqrt(h);
if (f > 0)
g = -g;
eivali[i] = scale * g;
eivali.coeffRef(i) = scale * g;
h = h - f * g;
eivalr[i-1] = f - g;
eivalr.coeffRef(i-1) = f - g;
eivali.start(i).setZero();
// Apply similarity transformation to remaining columns.
for (int j = 0; j < i; j++)
{
f = eivalr[j];
m_eivec(j,i) = f;
g = eivali[j] + m_eivec(j,j) * f;
f = eivalr.coeff(j);
m_eivec.coeffRef(j,i) = f;
g = eivali.coeff(j) + m_eivec.coeff(j,j) * f;
int bSize = i-j-1;
if (bSize>0)
{
g += (m_eivec.col(j).block(j+1, bSize).transpose() * eivalr.block(j+1, bSize))(0,0);
eivali.block(j+1, bSize) += m_eivec.col(j).block(j+1, bSize) * f;
}
eivali[j] = g;
eivali.coeffRef(j) = g;
}
f = (eivali.start(i).transpose() * eivalr.start(i))(0,0);
@ -206,15 +206,15 @@ void EigenSolver<MatrixType,IsSelfadjoint>::tridiagonalization(RealVectorType& e
eivalr.start(i) = m_eivec.row(i-1).start(i);
m_eivec.row(i).start(i).setZero();
}
eivalr[i] = h;
eivalr.coeffRef(i) = h;
}
// Accumulate transformations.
for (int i = 0; i < n-1; i++)
{
m_eivec(n-1,i) = m_eivec(i,i);
m_eivec(i,i) = 1.0;
Scalar h = eivalr[i+1];
m_eivec.coeffRef(n-1,i) = m_eivec.coeff(i,i);
m_eivec.coeffRef(i,i) = 1.0;
Scalar h = eivalr.coeff(i+1);
// FIXME this does not looks very stable ;)
if (h != 0.0)
{
@ -226,8 +226,8 @@ void EigenSolver<MatrixType,IsSelfadjoint>::tridiagonalization(RealVectorType& e
}
eivalr = m_eivec.row(eivalr.size()-1);
m_eivec.row(eivalr.size()-1).setZero();
m_eivec(n-1,n-1) = 1.0;
eivali[0] = 0.0;
m_eivec.coeffRef(n-1,n-1) = 1.0;
eivali.coeffRef(0) = 0.0;
}
@ -243,9 +243,9 @@ void EigenSolver<MatrixType,IsSelfadjoint>::tql2(RealVectorType& eivalr, RealVec
int n = eivalr.size();
for (int i = 1; i < n; i++) {
eivali[i-1] = eivali[i];
eivali.coeffRef(i-1) = eivali.coeff(i);
}
eivali[n-1] = 0.0;
eivali.coeffRef(n-1) = 0.0;
Scalar f = 0.0;
Scalar tst1 = 0.0;
@ -253,13 +253,13 @@ void EigenSolver<MatrixType,IsSelfadjoint>::tql2(RealVectorType& eivalr, RealVec
for (int l = 0; l < n; l++)
{
// Find small subdiagonal element
tst1 = std::max(tst1,ei_abs(eivalr[l]) + ei_abs(eivali[l]));
tst1 = std::max(tst1,ei_abs(eivalr.coeff(l)) + ei_abs(eivali.coeff(l)));
int m = l;
while ( (m < n) && (ei_abs(eivali[m]) > eps*tst1) )
while ( (m < n) && (ei_abs(eivali.coeff(m)) > eps*tst1) )
m++;
// If m == l, eivalr[l] is an eigenvalue,
// If m == l, eivalr.coeff(l) is an eigenvalue,
// otherwise, iterate.
if (m > l)
{
@ -269,26 +269,26 @@ void EigenSolver<MatrixType,IsSelfadjoint>::tql2(RealVectorType& eivalr, RealVec
iter = iter + 1;
// Compute implicit shift
Scalar g = eivalr[l];
Scalar p = (eivalr[l+1] - g) / (2.0 * eivali[l]);
Scalar g = eivalr.coeff(l);
Scalar p = (eivalr.coeff(l+1) - g) / (2.0 * eivali.coeff(l));
Scalar r = hypot(p,1.0);
if (p < 0)
r = -r;
eivalr[l] = eivali[l] / (p + r);
eivalr[l+1] = eivali[l] * (p + r);
Scalar dl1 = eivalr[l+1];
Scalar h = g - eivalr[l];
eivalr.coeffRef(l) = eivali.coeff(l) / (p + r);
eivalr.coeffRef(l+1) = eivali.coeff(l) * (p + r);
Scalar dl1 = eivalr.coeff(l+1);
Scalar h = g - eivalr.coeff(l);
if (l+2<n)
eivalr.end(n-l-2) -= RealVectorTypeX::constant(n-l-2, h);
f = f + h;
// Implicit QL transformation.
p = eivalr[m];
p = eivalr.coeff(m);
Scalar c = 1.0;
Scalar c2 = c;
Scalar c3 = c;
Scalar el1 = eivali[l+1];
Scalar el1 = eivali.coeff(l+1);
Scalar s = 0.0;
Scalar s2 = 0.0;
for (int i = m-1; i >= l; i--)
@ -296,32 +296,32 @@ void EigenSolver<MatrixType,IsSelfadjoint>::tql2(RealVectorType& eivalr, RealVec
c3 = c2;
c2 = c;
s2 = s;
g = c * eivali[i];
g = c * eivali.coeff(i);
h = c * p;
r = hypot(p,eivali[i]);
eivali[i+1] = s * r;
s = eivali[i] / r;
r = hypot(p,eivali.coeff(i));
eivali.coeffRef(i+1) = s * r;
s = eivali.coeff(i) / r;
c = p / r;
p = c * eivalr[i] - s * g;
eivalr[i+1] = h + s * (c * g + s * eivalr[i]);
p = c * eivalr.coeff(i) - s * g;
eivalr.coeffRef(i+1) = h + s * (c * g + s * eivalr.coeff(i));
// Accumulate transformation.
for (int k = 0; k < n; k++)
{
h = m_eivec(k,i+1);
m_eivec(k,i+1) = s * m_eivec(k,i) + c * h;
m_eivec(k,i) = c * m_eivec(k,i) - s * h;
h = m_eivec.coeff(k,i+1);
m_eivec.coeffRef(k,i+1) = s * m_eivec.coeff(k,i) + c * h;
m_eivec.coeffRef(k,i) = c * m_eivec.coeff(k,i) - s * h;
}
}
p = -s * s2 * c3 * el1 * eivali[l] / dl1;
eivali[l] = s * p;
eivalr[l] = c * p;
p = -s * s2 * c3 * el1 * eivali.coeff(l) / dl1;
eivali.coeffRef(l) = s * p;
eivalr.coeffRef(l) = c * p;
// Check for convergence.
} while (ei_abs(eivali[l]) > eps*tst1);
} while (ei_abs(eivali.coeff(l)) > eps*tst1);
}
eivalr[l] = eivalr[l] + f;
eivali[l] = 0.0;
eivalr.coeffRef(l) = eivalr.coeff(l) + f;
eivali.coeffRef(l) = 0.0;
}
// Sort eigenvalues and corresponding vectors.
@ -329,18 +329,18 @@ void EigenSolver<MatrixType,IsSelfadjoint>::tql2(RealVectorType& eivalr, RealVec
for (int i = 0; i < n-1; i++)
{
int k = i;
Scalar minValue = eivalr[i];
Scalar minValue = eivalr.coeff(i);
for (int j = i+1; j < n; j++)
{
if (eivalr[j] < minValue)
if (eivalr.coeff(j) < minValue)
{
k = j;
minValue = eivalr[j];
minValue = eivalr.coeff(j);
}
}
if (k != i)
{
std::swap(eivalr[i], eivalr[k]);
std::swap(eivalr.coeffRef(i), eivalr.coeffRef(k));
m_eivec.col(i).swap(m_eivec.col(k));
}
}
@ -371,14 +371,14 @@ void EigenSolver<MatrixType,IsSelfadjoint>::orthes(MatrixType& matH, RealVectorT
// FIXME could be rewritten, but this one looks better wrt cache
for (int i = high; i >= m; i--)
{
ort[i] = matH(i,m-1)/scale;
h += ort[i] * ort[i];
ort.coeffRef(i) = matH.coeff(i,m-1)/scale;
h += ort.coeff(i) * ort.coeff(i);
}
Scalar g = ei_sqrt(h);
if (ort[m] > 0)
if (ort.coeff(m) > 0)
g = -g;
h = h - ort[m] * g;
ort[m] = ort[m] - g;
h = h - ort.coeff(m) * g;
ort.coeffRef(m) = ort.coeff(m) - g;
// Apply Householder similarity transformation
// H = (I-u*u'/h)*H*(I-u*u')/h)
@ -389,8 +389,8 @@ void EigenSolver<MatrixType,IsSelfadjoint>::orthes(MatrixType& matH, RealVectorT
matH.block(0, m, high+1, bSize) -= ((matH.block(0, m, high+1, bSize) * ort.block(m, bSize)).lazy()
* (ort.block(m, bSize)/h).transpose()).lazy();
ort[m] = scale*ort[m];
matH(m,m-1) = scale*g;
ort.coeffRef(m) = scale*ort.coeff(m);
matH.coeffRef(m,m-1) = scale*g;
}
}
@ -399,12 +399,12 @@ void EigenSolver<MatrixType,IsSelfadjoint>::orthes(MatrixType& matH, RealVectorT
for (int m = high-1; m >= low+1; m--)
{
if (matH(m,m-1) != 0.0)
if (matH.coeff(m,m-1) != 0.0)
{
ort.block(m+1, high-m) = matH.col(m-1).block(m+1, high-m);
int bSize = high-m+1;
m_eivec.block(m, m, bSize, bSize) += ( (ort.block(m, bSize) / (matH(m,m-1) * ort[m] ) )
m_eivec.block(m, m, bSize, bSize) += ( (ort.block(m, bSize) / (matH.coeff(m,m-1) * ort.coeff(m) ) )
* (ort.block(m, bSize).transpose() * m_eivec.block(m, m, bSize, bSize)).lazy());
}
}
@ -458,8 +458,8 @@ void EigenSolver<MatrixType,IsSelfadjoint>::hqr2(MatrixType& matH)
// FIXME what's the purpose of the following since the condition is always false
if ((j < low) || (j > high))
{
m_eivalues[j].real() = matH(j,j);
m_eivalues[j].imag() = 0.0;
m_eivalues.coeffRef(j).real() = matH.coeff(j,j);
m_eivalues.coeffRef(j).imag() = 0.0;
}
norm += matH.col(j).start(std::min(j+1,nn)).cwiseAbs().sum();
}
@ -472,10 +472,10 @@ void EigenSolver<MatrixType,IsSelfadjoint>::hqr2(MatrixType& matH)
int l = n;
while (l > low)
{
s = ei_abs(matH(l-1,l-1)) + ei_abs(matH(l,l));
s = ei_abs(matH.coeff(l-1,l-1)) + ei_abs(matH.coeff(l,l));
if (s == 0.0)
s = norm;
if (ei_abs(matH(l,l-1)) < eps * s)
if (ei_abs(matH.coeff(l,l-1)) < eps * s)
break;
l--;
}
@ -484,21 +484,21 @@ void EigenSolver<MatrixType,IsSelfadjoint>::hqr2(MatrixType& matH)
// One root found
if (l == n)
{
matH(n,n) = matH(n,n) + exshift;
m_eivalues[n].real() = matH(n,n);
m_eivalues[n].imag() = 0.0;
matH.coeffRef(n,n) = matH.coeff(n,n) + exshift;
m_eivalues.coeffRef(n).real() = matH.coeff(n,n);
m_eivalues.coeffRef(n).imag() = 0.0;
n--;
iter = 0;
}
else if (l == n-1) // Two roots found
{
w = matH(n,n-1) * matH(n-1,n);
p = (matH(n-1,n-1) - matH(n,n)) / 2.0;
w = matH.coeff(n,n-1) * matH.coeff(n-1,n);
p = (matH.coeff(n-1,n-1) - matH.coeff(n,n)) / 2.0;
q = p * p + w;
z = ei_sqrt(ei_abs(q));
matH(n,n) = matH(n,n) + exshift;
matH(n-1,n-1) = matH(n-1,n-1) + exshift;
x = matH(n,n);
matH.coeffRef(n,n) = matH.coeff(n,n) + exshift;
matH.coeffRef(n-1,n-1) = matH.coeff(n-1,n-1) + exshift;
x = matH.coeff(n,n);
// Scalar pair
if (q >= 0)
@ -508,14 +508,14 @@ void EigenSolver<MatrixType,IsSelfadjoint>::hqr2(MatrixType& matH)
else
z = p - z;
m_eivalues[n-1].real() = x + z;
m_eivalues[n].real() = m_eivalues[n-1].real();
m_eivalues.coeffRef(n-1).real() = x + z;
m_eivalues.coeffRef(n).real() = m_eivalues.coeff(n-1).real();
if (z != 0.0)
m_eivalues[n].real() = x - w / z;
m_eivalues.coeffRef(n).real() = x - w / z;
m_eivalues[n-1].imag() = 0.0;
m_eivalues[n].imag() = 0.0;
x = matH(n,n-1);
m_eivalues.coeffRef(n-1).imag() = 0.0;
m_eivalues.coeffRef(n).imag() = 0.0;
x = matH.coeff(n,n-1);
s = ei_abs(x) + ei_abs(z);
p = x / s;
q = z / s;
@ -526,33 +526,33 @@ void EigenSolver<MatrixType,IsSelfadjoint>::hqr2(MatrixType& matH)
// Row modification
for (int j = n-1; j < nn; j++)
{
z = matH(n-1,j);
matH(n-1,j) = q * z + p * matH(n,j);
matH(n,j) = q * matH(n,j) - p * z;
z = matH.coeff(n-1,j);
matH.coeffRef(n-1,j) = q * z + p * matH.coeff(n,j);
matH.coeffRef(n,j) = q * matH.coeff(n,j) - p * z;
}
// Column modification
for (int i = 0; i <= n; i++)
{
z = matH(i,n-1);
matH(i,n-1) = q * z + p * matH(i,n);
matH(i,n) = q * matH(i,n) - p * z;
z = matH.coeff(i,n-1);
matH.coeffRef(i,n-1) = q * z + p * matH.coeff(i,n);
matH.coeffRef(i,n) = q * matH.coeff(i,n) - p * z;
}
// Accumulate transformations
for (int i = low; i <= high; i++)
{
z = m_eivec(i,n-1);
m_eivec(i,n-1) = q * z + p * m_eivec(i,n);
m_eivec(i,n) = q * m_eivec(i,n) - p * z;
z = m_eivec.coeff(i,n-1);
m_eivec.coeffRef(i,n-1) = q * z + p * m_eivec.coeff(i,n);
m_eivec.coeffRef(i,n) = q * m_eivec.coeff(i,n) - p * z;
}
}
else // Complex pair
{
m_eivalues[n-1].real() = x + p;
m_eivalues[n].real() = x + p;
m_eivalues[n-1].imag() = z;
m_eivalues[n].imag() = -z;
m_eivalues.coeffRef(n-1).real() = x + p;
m_eivalues.coeffRef(n).real() = x + p;
m_eivalues.coeffRef(n-1).imag() = z;
m_eivalues.coeffRef(n).imag() = -z;
}
n = n - 2;
iter = 0;
@ -560,13 +560,13 @@ void EigenSolver<MatrixType,IsSelfadjoint>::hqr2(MatrixType& matH)
else // No convergence yet
{
// Form shift
x = matH(n,n);
x = matH.coeff(n,n);
y = 0.0;
w = 0.0;
if (l < n)
{
y = matH(n-1,n-1);
w = matH(n,n-1) * matH(n-1,n);
y = matH.coeff(n-1,n-1);
w = matH.coeff(n,n-1) * matH.coeff(n-1,n);
}
// Wilkinson's original ad hoc shift
@ -574,8 +574,8 @@ void EigenSolver<MatrixType,IsSelfadjoint>::hqr2(MatrixType& matH)
{
exshift += x;
for (int i = low; i <= n; i++)
matH(i,i) -= x;
s = ei_abs(matH(n,n-1)) + ei_abs(matH(n-1,n-2));
matH.coeffRef(i,i) -= x;
s = ei_abs(matH.coeff(n,n-1)) + ei_abs(matH.coeff(n-1,n-2));
x = y = 0.75 * s;
w = -0.4375 * s * s;
}
@ -592,7 +592,7 @@ void EigenSolver<MatrixType,IsSelfadjoint>::hqr2(MatrixType& matH)
s = -s;
s = x - w / ((y - x) / 2.0 + s);
for (int i = low; i <= n; i++)
matH(i,i) -= s;
matH.coeffRef(i,i) -= s;
exshift += s;
x = y = w = 0.964;
}
@ -604,12 +604,12 @@ void EigenSolver<MatrixType,IsSelfadjoint>::hqr2(MatrixType& matH)
int m = n-2;
while (m >= l)
{
z = matH(m,m);
z = matH.coeff(m,m);
r = x - z;
s = y - z;
p = (r * s - w) / matH(m+1,m) + matH(m,m+1);
q = matH(m+1,m+1) - z - r - s;
r = matH(m+2,m+1);
p = (r * s - w) / matH.coeff(m+1,m) + matH.coeff(m,m+1);
q = matH.coeff(m+1,m+1) - z - r - s;
r = matH.coeff(m+2,m+1);
s = ei_abs(p) + ei_abs(q) + ei_abs(r);
p = p / s;
q = q / s;
@ -617,9 +617,9 @@ void EigenSolver<MatrixType,IsSelfadjoint>::hqr2(MatrixType& matH)
if (m == l) {
break;
}
if (ei_abs(matH(m,m-1)) * (ei_abs(q) + ei_abs(r)) <
eps * (ei_abs(p) * (ei_abs(matH(m-1,m-1)) + ei_abs(z) +
ei_abs(matH(m+1,m+1)))))
if (ei_abs(matH.coeff(m,m-1)) * (ei_abs(q) + ei_abs(r)) <
eps * (ei_abs(p) * (ei_abs(matH.coeff(m-1,m-1)) + ei_abs(z) +
ei_abs(matH.coeff(m+1,m+1)))))
{
break;
}
@ -628,9 +628,9 @@ void EigenSolver<MatrixType,IsSelfadjoint>::hqr2(MatrixType& matH)
for (int i = m+2; i <= n; i++)
{
matH(i,i-2) = 0.0;
matH.coeffRef(i,i-2) = 0.0;
if (i > m+2)
matH(i,i-3) = 0.0;
matH.coeffRef(i,i-3) = 0.0;
}
// Double QR step involving rows l:n and columns m:n
@ -638,9 +638,9 @@ void EigenSolver<MatrixType,IsSelfadjoint>::hqr2(MatrixType& matH)
{
int notlast = (k != n-1);
if (k != m) {
p = matH(k,k-1);
q = matH(k+1,k-1);
r = (notlast ? matH(k+2,k-1) : 0.0);
p = matH.coeff(k,k-1);
q = matH.coeff(k+1,k-1);
r = (notlast ? matH.coeff(k+2,k-1) : 0.0);
x = ei_abs(p) + ei_abs(q) + ei_abs(r);
if (x != 0.0)
{
@ -661,9 +661,9 @@ void EigenSolver<MatrixType,IsSelfadjoint>::hqr2(MatrixType& matH)
if (s != 0)
{
if (k != m)
matH(k,k-1) = -s * x;
matH.coeffRef(k,k-1) = -s * x;
else if (l != m)
matH(k,k-1) = -matH(k,k-1);
matH.coeffRef(k,k-1) = -matH.coeff(k,k-1);
p = p + s;
x = p / s;
@ -675,40 +675,40 @@ void EigenSolver<MatrixType,IsSelfadjoint>::hqr2(MatrixType& matH)
// Row modification
for (int j = k; j < nn; j++)
{
p = matH(k,j) + q * matH(k+1,j);
p = matH.coeff(k,j) + q * matH.coeff(k+1,j);
if (notlast)
{
p = p + r * matH(k+2,j);
matH(k+2,j) = matH(k+2,j) - p * z;
p = p + r * matH.coeff(k+2,j);
matH.coeffRef(k+2,j) = matH.coeff(k+2,j) - p * z;
}
matH(k,j) = matH(k,j) - p * x;
matH(k+1,j) = matH(k+1,j) - p * y;
matH.coeffRef(k,j) = matH.coeff(k,j) - p * x;
matH.coeffRef(k+1,j) = matH.coeff(k+1,j) - p * y;
}
// Column modification
for (int i = 0; i <= std::min(n,k+3); i++)
{
p = x * matH(i,k) + y * matH(i,k+1);
p = x * matH.coeff(i,k) + y * matH.coeff(i,k+1);
if (notlast)
{
p = p + z * matH(i,k+2);
matH(i,k+2) = matH(i,k+2) - p * r;
p = p + z * matH.coeff(i,k+2);
matH.coeffRef(i,k+2) = matH.coeff(i,k+2) - p * r;
}
matH(i,k) = matH(i,k) - p;
matH(i,k+1) = matH(i,k+1) - p * q;
matH.coeffRef(i,k) = matH.coeff(i,k) - p;
matH.coeffRef(i,k+1) = matH.coeff(i,k+1) - p * q;
}
// Accumulate transformations
for (int i = low; i <= high; i++)
{
p = x * m_eivec(i,k) + y * m_eivec(i,k+1);
p = x * m_eivec.coeff(i,k) + y * m_eivec.coeff(i,k+1);
if (notlast)
{
p = p + z * m_eivec(i,k+2);
m_eivec(i,k+2) = m_eivec(i,k+2) - p * r;
p = p + z * m_eivec.coeff(i,k+2);
m_eivec.coeffRef(i,k+2) = m_eivec.coeff(i,k+2) - p * r;
}
m_eivec(i,k) = m_eivec(i,k) - p;
m_eivec(i,k+1) = m_eivec(i,k+1) - p * q;
m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) - p;
m_eivec.coeffRef(i,k+1) = m_eivec.coeff(i,k+1) - p * q;
}
} // (s != 0)
} // k loop
@ -723,20 +723,20 @@ void EigenSolver<MatrixType,IsSelfadjoint>::hqr2(MatrixType& matH)
for (n = nn-1; n >= 0; n--)
{
p = m_eivalues[n].real();
q = m_eivalues[n].imag();
p = m_eivalues.coeff(n).real();
q = m_eivalues.coeff(n).imag();
// Scalar vector
if (q == 0)
{
int l = n;
matH(n,n) = 1.0;
matH.coeffRef(n,n) = 1.0;
for (int i = n-1; i >= 0; i--)
{
w = matH(i,i) - p;
w = matH.coeff(i,i) - p;
r = (matH.row(i).end(nn-l) * matH.col(n).end(nn-l))(0,0);
if (m_eivalues[i].imag() < 0.0)
if (m_eivalues.coeff(i).imag() < 0.0)
{
z = w;
s = r;
@ -744,28 +744,28 @@ void EigenSolver<MatrixType,IsSelfadjoint>::hqr2(MatrixType& matH)
else
{
l = i;
if (m_eivalues[i].imag() == 0.0)
if (m_eivalues.coeff(i).imag() == 0.0)
{
if (w != 0.0)
matH(i,n) = -r / w;
matH.coeffRef(i,n) = -r / w;
else
matH(i,n) = -r / (eps * norm);
matH.coeffRef(i,n) = -r / (eps * norm);
}
else // Solve real equations
{
x = matH(i,i+1);
y = matH(i+1,i);
q = (m_eivalues[i].real() - p) * (m_eivalues[i].real() - p) + m_eivalues[i].imag() * m_eivalues[i].imag();
x = matH.coeff(i,i+1);
y = matH.coeff(i+1,i);
q = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
t = (x * s - z * r) / q;
matH(i,n) = t;
matH.coeffRef(i,n) = t;
if (ei_abs(x) > ei_abs(z))
matH(i+1,n) = (-r - w * t) / x;
matH.coeffRef(i+1,n) = (-r - w * t) / x;
else
matH(i+1,n) = (-s - y * t) / z;
matH.coeffRef(i+1,n) = (-s - y * t) / z;
}
// Overflow control
t = ei_abs(matH(i,n));
t = ei_abs(matH.coeff(i,n));
if ((eps * t) * t > 1)
matH.col(n).end(nn-i) /= t;
}
@ -777,27 +777,27 @@ void EigenSolver<MatrixType,IsSelfadjoint>::hqr2(MatrixType& matH)
int l = n-1;
// Last vector component imaginary so matrix is triangular
if (ei_abs(matH(n,n-1)) > ei_abs(matH(n-1,n)))
if (ei_abs(matH.coeff(n,n-1)) > ei_abs(matH.coeff(n-1,n)))
{
matH(n-1,n-1) = q / matH(n,n-1);
matH(n-1,n) = -(matH(n,n) - p) / matH(n,n-1);
matH.coeffRef(n-1,n-1) = q / matH.coeff(n,n-1);
matH.coeffRef(n-1,n) = -(matH.coeff(n,n) - p) / matH.coeff(n,n-1);
}
else
{
cc = cdiv<Scalar>(0.0,-matH(n-1,n),matH(n-1,n-1)-p,q);
matH(n-1,n-1) = ei_real(cc);
matH(n-1,n) = ei_imag(cc);
cc = cdiv<Scalar>(0.0,-matH.coeff(n-1,n),matH.coeff(n-1,n-1)-p,q);
matH.coeffRef(n-1,n-1) = ei_real(cc);
matH.coeffRef(n-1,n) = ei_imag(cc);
}
matH(n,n-1) = 0.0;
matH(n,n) = 1.0;
matH.coeffRef(n,n-1) = 0.0;
matH.coeffRef(n,n) = 1.0;
for (int i = n-2; i >= 0; i--)
{
Scalar ra,sa,vr,vi;
ra = (matH.row(i).end(nn-l) * matH.col(n-1).end(nn-l)).lazy()(0,0);
sa = (matH.row(i).end(nn-l) * matH.col(n).end(nn-l)).lazy()(0,0);
w = matH(i,i) - p;
w = matH.coeff(i,i) - p;
if (m_eivalues[i].imag() < 0.0)
if (m_eivalues.coeff(i).imag() < 0.0)
{
z = w;
r = ra;
@ -806,40 +806,40 @@ void EigenSolver<MatrixType,IsSelfadjoint>::hqr2(MatrixType& matH)
else
{
l = i;
if (m_eivalues[i].imag() == 0)
if (m_eivalues.coeff(i).imag() == 0)
{
cc = cdiv(-ra,-sa,w,q);
matH(i,n-1) = ei_real(cc);
matH(i,n) = ei_imag(cc);
matH.coeffRef(i,n-1) = ei_real(cc);
matH.coeffRef(i,n) = ei_imag(cc);
}
else
{
// Solve complex equations
x = matH(i,i+1);
y = matH(i+1,i);
vr = (m_eivalues[i].real() - p) * (m_eivalues[i].real() - p) + m_eivalues[i].imag() * m_eivalues[i].imag() - q * q;
vi = (m_eivalues[i].real() - p) * 2.0 * q;
x = matH.coeff(i,i+1);
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;
vi = (m_eivalues.coeff(i).real() - p) * 2.0 * q;
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));
cc= cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi);
matH(i,n-1) = ei_real(cc);
matH(i,n) = ei_imag(cc);
matH.coeffRef(i,n-1) = ei_real(cc);
matH.coeffRef(i,n) = ei_imag(cc);
if (ei_abs(x) > (ei_abs(z) + ei_abs(q)))
{
matH(i+1,n-1) = (-ra - w * matH(i,n-1) + q * matH(i,n)) / x;
matH(i+1,n) = (-sa - w * matH(i,n) - q * matH(i,n-1)) / x;
matH.coeffRef(i+1,n-1) = (-ra - w * matH.coeff(i,n-1) + q * matH.coeff(i,n)) / x;
matH.coeffRef(i+1,n) = (-sa - w * matH.coeff(i,n) - q * matH.coeff(i,n-1)) / x;
}
else
{
cc = cdiv(-r-y*matH(i,n-1),-s-y*matH(i,n),z,q);
matH(i+1,n-1) = ei_real(cc);
matH(i+1,n) = ei_imag(cc);
cc = cdiv(-r-y*matH.coeff(i,n-1),-s-y*matH.coeff(i,n),z,q);
matH.coeffRef(i+1,n-1) = ei_real(cc);
matH.coeffRef(i+1,n) = ei_imag(cc);
}
}
// Overflow control
t = std::max(ei_abs(matH(i,n-1)),ei_abs(matH(i,n)));
t = std::max(ei_abs(matH.coeff(i,n-1)),ei_abs(matH.coeff(i,n)));
if ((eps * t) * t > 1)
matH.block(i, n-1, nn-i, 2) /= t;