From d9c1224133153045fb92df351bdd3660a801f8d7 Mon Sep 17 00:00:00 2001 From: Jitse Niesen Date: Mon, 26 Apr 2010 17:43:16 +0100 Subject: Simplify computation of eigenvectors in EigenSolver. * reduce scope of declarations * use that low = 0 and high = size-1 * rename some variables * rename hqr2_step2() to computeEigenvectors() * exploit that ei_isMuchSmallerThan takes absolute value of arguments --- Eigen/src/Eigenvalues/EigenSolver.h | 107 ++++++++++++++++-------------------- 1 file changed, 46 insertions(+), 61 deletions(-) diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h index 114bfab6f..f8b953d9b 100644 --- a/Eigen/src/Eigenvalues/EigenSolver.h +++ b/Eigen/src/Eigenvalues/EigenSolver.h @@ -261,7 +261,7 @@ template class EigenSolver EigenSolver& compute(const MatrixType& matrix); private: - void hqr2_step2(MatrixType& matH); + void computeEigenvectors(MatrixType& matH); protected: MatrixType m_eivec; @@ -297,7 +297,7 @@ typename EigenSolver::EigenvectorsType EigenSolver::eige EigenvectorsType matV(n,n); for (int j=0; j(); @@ -349,7 +349,7 @@ EigenSolver& EigenSolver::compute(const MatrixType& matr } // Compute eigenvectors. - hqr2_step2(matT); + computeEigenvectors(matT); m_isInitialized = true; return *this; @@ -376,19 +376,16 @@ std::complex cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi) template -void EigenSolver::hqr2_step2(MatrixType& matH) +void EigenSolver::computeEigenvectors(MatrixType& matH) { - const int nn = m_eivec.cols(); - const int low = 0; - const int high = nn-1; - const Scalar eps = ei_pow(Scalar(2),ei_is_same_type::ret ? Scalar(-23) : Scalar(-52)); - Scalar p, q, r=0, s=0, t, w, x, y, z=0; + const int size = m_eivec.cols(); + const Scalar eps = NumTraits::epsilon(); // inefficient! this is already computed in RealSchur Scalar norm = 0.0; - for (int j = 0; j < nn; ++j) + for (int j = 0; j < size; ++j) { - norm += matH.row(j).segment(std::max(j-1,0), nn-std::max(j-1,0)).cwiseAbs().sum(); + norm += matH.row(j).segment(std::max(j-1,0), size-std::max(j-1,0)).cwiseAbs().sum(); } // Backsubstitute to find vectors of upper triangular form @@ -397,25 +394,27 @@ void EigenSolver::hqr2_step2(MatrixType& matH) return; } - for (int n = nn-1; n >= 0; n--) + for (int n = size-1; n >= 0; n--) { - p = m_eivalues.coeff(n).real(); - q = m_eivalues.coeff(n).imag(); + Scalar p = m_eivalues.coeff(n).real(); + Scalar q = m_eivalues.coeff(n).imag(); // Scalar vector if (q == 0) { + Scalar lastr=0, lastw=0; int l = n; + matH.coeffRef(n,n) = 1.0; for (int i = n-1; i >= 0; i--) { - w = matH.coeff(i,i) - p; - r = matH.row(i).segment(l,n-l+1).dot(matH.col(n).segment(l, n-l+1)); + Scalar w = matH.coeff(i,i) - p; + Scalar r = matH.row(i).segment(l,n-l+1).dot(matH.col(n).segment(l, n-l+1)); if (m_eivalues.coeff(i).imag() < 0.0) { - z = w; - s = r; + lastw = w; + lastr = r; } else { @@ -429,27 +428,27 @@ void EigenSolver::hqr2_step2(MatrixType& matH) } else // Solve real equations { - 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; + Scalar x = matH.coeff(i,i+1); + Scalar y = matH.coeff(i+1,i); + Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag(); + Scalar t = (x * lastr - lastw * r) / denom; matH.coeffRef(i,n) = t; - if (ei_abs(x) > ei_abs(z)) + if (ei_abs(x) > ei_abs(lastw)) matH.coeffRef(i+1,n) = (-r - w * t) / x; else - matH.coeffRef(i+1,n) = (-s - y * t) / z; + matH.coeffRef(i+1,n) = (-lastr - y * t) / lastw; } // Overflow control - t = ei_abs(matH.coeff(i,n)); + Scalar t = ei_abs(matH.coeff(i,n)); if ((eps * t) * t > 1) - matH.col(n).tail(nn-i) /= t; + matH.col(n).tail(size-i) /= t; } } } else if (q < 0) // Complex vector { - std::complex cc; + Scalar lastra=0, lastsa=0, lastw=0; int l = n-1; // Last vector component imaginary so matrix is triangular @@ -460,7 +459,7 @@ void EigenSolver::hqr2_step2(MatrixType& matH) } else { - cc = cdiv(0.0,-matH.coeff(n-1,n),matH.coeff(n-1,n-1)-p,q); + std::complex cc = cdiv(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); } @@ -468,79 +467,65 @@ void EigenSolver::hqr2_step2(MatrixType& matH) matH.coeffRef(n,n) = 1.0; for (int i = n-2; i >= 0; i--) { - Scalar ra,sa,vr,vi; - ra = matH.row(i).segment(l, n-l+1).dot(matH.col(n-1).segment(l, n-l+1)); - sa = matH.row(i).segment(l, n-l+1).dot(matH.col(n).segment(l, n-l+1)); - w = matH.coeff(i,i) - p; + Scalar ra = matH.row(i).segment(l, n-l+1).dot(matH.col(n-1).segment(l, n-l+1)); + Scalar sa = matH.row(i).segment(l, n-l+1).dot(matH.col(n).segment(l, n-l+1)); + Scalar w = matH.coeff(i,i) - p; if (m_eivalues.coeff(i).imag() < 0.0) { - z = w; - r = ra; - s = sa; + lastw = w; + lastra = ra; + lastsa = sa; } else { l = i; if (m_eivalues.coeff(i).imag() == 0) { - cc = cdiv(-ra,-sa,w,q); + std::complex cc = cdiv(-ra,-sa,w,q); matH.coeffRef(i,n-1) = ei_real(cc); matH.coeffRef(i,n) = ei_imag(cc); } else { // Solve complex equations - 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) * Scalar(2) * q; + Scalar x = matH.coeff(i,i+1); + Scalar y = matH.coeff(i+1,i); + Scalar 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; + Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * 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)); + vr = eps * norm * (ei_abs(w) + ei_abs(q) + ei_abs(x) + ei_abs(y) + ei_abs(lastw)); - cc= cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi); + std::complex cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi); 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))) + if (ei_abs(x) > (ei_abs(lastw) + ei_abs(q))) { 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.coeff(i,n-1),-s-y*matH.coeff(i,n),z,q); + cc = cdiv(-lastra-y*matH.coeff(i,n-1),-lastsa-y*matH.coeff(i,n),lastw,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.coeff(i,n-1)),ei_abs(matH.coeff(i,n))); + Scalar 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; + matH.block(i, n-1, size-i, 2) /= t; } } } } - // Vectors of isolated roots - for (int i = 0; i < nn; ++i) - { - // FIXME again what's the purpose of this test ? - // in this algo low==0 and high==nn-1 !! - if (i < low || i > high) - { - m_eivec.row(i).tail(nn-i) = matH.row(i).tail(nn-i); - } - } - // Back transformation to get eigenvectors of original matrix - int bRows = high-low+1; - for (int j = nn-1; j >= low; j--) + for (int j = size-1; j >= 0; j--) { - int bSize = std::min(j,high)-low+1; - m_eivec.col(j).segment(low, bRows) = (m_eivec.block(low, low, bRows, bSize) * matH.col(j).segment(low, bSize)); + m_eivec.col(j).segment(0, size) = m_eivec.leftCols(j+1) * matH.col(j).segment(0, j+1); } } -- cgit v1.2.3