diff options
author | Gael Guennebaud <g.gael@free.fr> | 2009-07-31 17:35:20 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2009-07-31 17:35:20 +0200 |
commit | 2796bcabb1151ad8de2bf2ab9117baea40ae4d30 (patch) | |
tree | 94bbdb7be79823960042da2d565cb8c5473e6359 /Eigen | |
parent | a156f5a8696a22a2738ab54bcc77a7565d5e9019 (diff) |
some cleaning
Diffstat (limited to 'Eigen')
-rw-r--r-- | Eigen/src/QR/SelfAdjointEigenSolver.h | 30 | ||||
-rw-r--r-- | Eigen/src/QR/Tridiagonalization.h | 10 |
2 files changed, 15 insertions, 25 deletions
diff --git a/Eigen/src/QR/SelfAdjointEigenSolver.h b/Eigen/src/QR/SelfAdjointEigenSolver.h index 167fedebf..b003fb4d7 100644 --- a/Eigen/src/QR/SelfAdjointEigenSolver.h +++ b/Eigen/src/QR/SelfAdjointEigenSolver.h @@ -383,24 +383,18 @@ static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int st int kn1 = (k+1)*n; #endif // let's do the product manually to avoid the need of temporaries... - Matrix<Scalar,Dynamic,1> aux = Map<Matrix<Scalar,Dynamic,1>, ForceAligned >(matrixQ+kn,n); - Map<Matrix<Scalar,Dynamic,1>, ForceAligned >(matrixQ+kn,n) - = Map<Matrix<Scalar,Dynamic,1>, ForceAligned >(matrixQ+kn,n) * c - s * Map<Matrix<Scalar,Dynamic,1> >(matrixQ+kn1,n); - - Map<Matrix<Scalar,Dynamic,1>, ForceAligned >(matrixQ+kn1,n) - = Map<Matrix<Scalar,Dynamic,1>, ForceAligned >(matrixQ+kn1,n) * c - s * aux;//Map<Matrix<Scalar,Dynamic,1> >(matrixQ+kn,n); -// for (int i=0; i<n; ++i) -// { -// #ifdef EIGEN_DEFAULT_TO_ROW_MAJOR -// Scalar matrixQ_i_k = matrixQ[i*n+k]; -// matrixQ[i*n+k] = c * matrixQ_i_k - s * matrixQ[i*n+k+1]; -// matrixQ[i*n+k+1] = s * matrixQ_i_k + c * matrixQ[i*n+k+1]; -// #else -// Scalar matrixQ_i_k = matrixQ[i+kn]; -// matrixQ[i+kn] = c * matrixQ_i_k - s * matrixQ[i+kn1]; -// matrixQ[i+kn1] = s * matrixQ_i_k + c * matrixQ[i+kn1]; -// #endif -// } + for (int i=0; i<n; ++i) + { + #ifdef EIGEN_DEFAULT_TO_ROW_MAJOR + Scalar matrixQ_i_k = matrixQ[i*n+k]; + matrixQ[i*n+k] = c * matrixQ_i_k - s * matrixQ[i*n+k+1]; + matrixQ[i*n+k+1] = s * matrixQ_i_k + c * matrixQ[i*n+k+1]; + #else + Scalar matrixQ_i_k = matrixQ[i+kn]; + matrixQ[i+kn] = c * matrixQ_i_k - s * matrixQ[i+kn1]; + matrixQ[i+kn1] = s * matrixQ_i_k + c * matrixQ[i+kn1]; + #endif + } } } } diff --git a/Eigen/src/QR/Tridiagonalization.h b/Eigen/src/QR/Tridiagonalization.h index d8d537ad0..e6f594237 100644 --- a/Eigen/src/QR/Tridiagonalization.h +++ b/Eigen/src/QR/Tridiagonalization.h @@ -226,12 +226,8 @@ void Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& // i.e., A = H' A H where H = I - h v v' and v = matA.col(i).end(n-i-1) matA.col(i).coeffRef(i+1) = 1; -// hCoeffs.end(n-i-1) = (matA.corner(BottomRight,n-i-1,n-i-1).template selfadjointView<LowerTriangular>() -// * (h * matA.col(i).end(n-i-1))); - - hCoeffs.end(n-i-1).setZero(); - ei_product_selfadjoint_vector<Scalar,MatrixType::Flags&RowMajorBit,LowerTriangular,false,false> - (n-i-1,matA.corner(BottomRight,n-i-1,n-i-1).data(), matA.stride(), matA.col(i).end(n-i-1).data(), 1, const_cast<Scalar*>(hCoeffs.end(n-i-1).data()), h); + hCoeffs.end(n-i-1) = (matA.corner(BottomRight,n-i-1,n-i-1).template selfadjointView<LowerTriangular>() + * (h * matA.col(i).end(n-i-1))); hCoeffs.end(n-i-1) += (h*Scalar(-0.5)*(matA.col(i).end(n-i-1).dot(hCoeffs.end(n-i-1)))) * matA.col(i).end(n-i-1); @@ -290,7 +286,7 @@ void Tridiagonalization<MatrixType>::matrixQInPlace(MatrixBase<QDerived>* q) con aux.end(n-i-1) = (m_hCoeffs.coeff(i) * m_matrix.col(i).end(n-i-1).adjoint() * matQ.corner(BottomRight,n-i-1,n-i-1)).lazy(); // rank one update, TODO ! make it works efficiently as expected for (int j=i+1;j<n;++j) - matQ.col(j).end(n-i-1) -= ( aux.coeff(j)) * m_matrix.col(i).end(n-i-1); + matQ.col(j).end(n-i-1) -= aux.coeff(j) * m_matrix.col(i).end(n-i-1); // matQ.corner(BottomRight,n-i-1,n-i-1) -= (m_matrix.col(i).end(n-i-1) * aux.end(n-i-1)).lazy(); m_matrix.const_cast_derived().coeffRef(i+1,i) = tmp; |