diff options
author | Gael Guennebaud <g.gael@free.fr> | 2008-05-31 16:31:10 +0000 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2008-05-31 16:31:10 +0000 |
commit | a2f71f9d7e1e443814fe60726d99a4b0508baefa (patch) | |
tree | 1d5ad7a09fcfcd32fc04b5b450a1fa5b19f4916b | |
parent | c9fb248c3667af2a9fbac3011a723f1ec32f1601 (diff) |
updated EigenSolver to use .coeff / .coeffRef
-rw-r--r-- | Eigen/src/QR/EigenSolver.h | 322 |
1 files changed, 161 insertions, 161 deletions
diff --git a/Eigen/src/QR/EigenSolver.h b/Eigen/src/QR/EigenSolver.h index cf3ea9c94..4e875757e 100644 --- a/Eigen/src/QR/EigenSolver.h +++ b/Eigen/src/QR/EigenSolver.h @@ -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; |