mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-12-15 07:10:37 +08:00
Code Issues Packages Projects Releases Wiki Activity

some cleaning

Browse Source
This commit is contained in:
Gael Guennebaud 2009-07-31 17:35:20 +02:00
parent a156f5a869
commit 2796bcabb1
2 changed files with 15 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
Eigen/src/QR/SelfAdjointEigenSolver.h
View File

@ -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
}
}
}
}

10
Eigen/src/QR/Tridiagonalization.h
View File

@ -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;