aboutsummaryrefslogtreecommitdiffhomepage
path: root/Eigen
diff options
context:
space:
mode:
authorGravatar Gael Guennebaud <g.gael@free.fr>2009-07-31 17:35:20 +0200
committerGravatar Gael Guennebaud <g.gael@free.fr>2009-07-31 17:35:20 +0200
commit2796bcabb1151ad8de2bf2ab9117baea40ae4d30 (patch)
tree94bbdb7be79823960042da2d565cb8c5473e6359 /Eigen
parenta156f5a8696a22a2738ab54bcc77a7565d5e9019 (diff)
some cleaning
Diffstat (limited to 'Eigen')
-rw-r--r--Eigen/src/QR/SelfAdjointEigenSolver.h30
-rw-r--r--Eigen/src/QR/Tridiagonalization.h10
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;