mirror of
https://gitlab.com/libeigen/eigen.git
synced 2026-04-10 11:34:33 +08:00
bug #86 : use internal:: namespace instead of ei_ prefix
This commit is contained in:
Benoit Jacob
@@ -169,8 +169,8 @@ template<typename _MatrixType> class ComplexEigenSolver
|
||||
*/
|
||||
const EigenvectorType& eigenvectors() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
|
||||
ei_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
|
||||
eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
|
||||
eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
|
||||
return m_eivec;
|
||||
}
|
||||
|
||||
@@ -193,7 +193,7 @@ template<typename _MatrixType> class ComplexEigenSolver
|
||||
*/
|
||||
const EigenvalueType& eigenvalues() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
|
||||
eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
|
||||
return m_eivalues;
|
||||
}
|
||||
|
||||
@@ -229,7 +229,7 @@ template<typename _MatrixType> class ComplexEigenSolver
|
||||
*/
|
||||
ComputationInfo info() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
|
||||
eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
|
||||
return m_schur.info();
|
||||
}
|
||||
|
||||
@@ -293,7 +293,7 @@ void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm
|
||||
{
|
||||
// If the i-th and k-th eigenvalue are equal, then z equals 0.
|
||||
// Use a small value instead, to prevent division by zero.
|
||||
ei_real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
|
||||
internal::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
|
||||
}
|
||||
m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
|
||||
}
|
||||
|
||||
@@ -30,7 +30,9 @@
|
||||
#include "./EigenvaluesCommon.h"
|
||||
#include "./HessenbergDecomposition.h"
|
||||
|
||||
template<typename MatrixType, bool IsComplex> struct ei_complex_schur_reduce_to_hessenberg;
|
||||
namespace internal {
|
||||
template<typename MatrixType, bool IsComplex> struct complex_schur_reduce_to_hessenberg;
|
||||
}
|
||||
|
||||
/** \eigenvalues_module \ingroup Eigenvalues_Module
|
||||
*
|
||||
@@ -146,8 +148,8 @@ template<typename _MatrixType> class ComplexSchur
|
||||
*/
|
||||
const ComplexMatrixType& matrixU() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "ComplexSchur is not initialized.");
|
||||
ei_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition.");
|
||||
eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
|
||||
eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition.");
|
||||
return m_matU;
|
||||
}
|
||||
|
||||
@@ -170,7 +172,7 @@ template<typename _MatrixType> class ComplexSchur
|
||||
*/
|
||||
const ComplexMatrixType& matrixT() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "ComplexSchur is not initialized.");
|
||||
eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
|
||||
return m_matT;
|
||||
}
|
||||
|
||||
@@ -201,7 +203,7 @@ template<typename _MatrixType> class ComplexSchur
|
||||
*/
|
||||
ComputationInfo info() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "RealSchur is not initialized.");
|
||||
eigen_assert(m_isInitialized && "RealSchur is not initialized.");
|
||||
return m_info;
|
||||
}
|
||||
|
||||
@@ -222,22 +224,24 @@ template<typename _MatrixType> class ComplexSchur
|
||||
bool subdiagonalEntryIsNeglegible(Index i);
|
||||
ComplexScalar computeShift(Index iu, Index iter);
|
||||
void reduceToTriangularForm(bool computeU);
|
||||
friend struct ei_complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>;
|
||||
friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>;
|
||||
};
|
||||
|
||||
namespace internal {
|
||||
|
||||
/** Computes the principal value of the square root of the complex \a z. */
|
||||
template<typename RealScalar>
|
||||
std::complex<RealScalar> ei_sqrt(const std::complex<RealScalar> &z)
|
||||
std::complex<RealScalar> sqrt(const std::complex<RealScalar> &z)
|
||||
{
|
||||
RealScalar t, tre, tim;
|
||||
|
||||
t = ei_abs(z);
|
||||
t = abs(z);
|
||||
|
||||
if (ei_abs(ei_real(z)) <= ei_abs(ei_imag(z)))
|
||||
if (abs(real(z)) <= abs(imag(z)))
|
||||
{
|
||||
// No cancellation in these formulas
|
||||
tre = ei_sqrt(RealScalar(0.5)*(t + ei_real(z)));
|
||||
tim = ei_sqrt(RealScalar(0.5)*(t - ei_real(z)));
|
||||
tre = sqrt(RealScalar(0.5)*(t + real(z)));
|
||||
tim = sqrt(RealScalar(0.5)*(t - real(z)));
|
||||
}
|
||||
else
|
||||
{
|
||||
@@ -245,14 +249,14 @@ std::complex<RealScalar> ei_sqrt(const std::complex<RealScalar> &z)
|
||||
if (z.real() > RealScalar(0))
|
||||
{
|
||||
tre = t + z.real();
|
||||
tim = ei_abs(ei_imag(z))*ei_sqrt(RealScalar(0.5)/tre);
|
||||
tre = ei_sqrt(RealScalar(0.5)*tre);
|
||||
tim = abs(imag(z))*sqrt(RealScalar(0.5)/tre);
|
||||
tre = sqrt(RealScalar(0.5)*tre);
|
||||
}
|
||||
else
|
||||
{
|
||||
tim = t - z.real();
|
||||
tre = ei_abs(ei_imag(z))*ei_sqrt(RealScalar(0.5)/tim);
|
||||
tim = ei_sqrt(RealScalar(0.5)*tim);
|
||||
tre = abs(imag(z))*sqrt(RealScalar(0.5)/tim);
|
||||
tim = sqrt(RealScalar(0.5)*tim);
|
||||
}
|
||||
}
|
||||
if(z.imag() < RealScalar(0))
|
||||
@@ -260,6 +264,7 @@ std::complex<RealScalar> ei_sqrt(const std::complex<RealScalar> &z)
|
||||
|
||||
return (std::complex<RealScalar>(tre,tim));
|
||||
}
|
||||
} // end namespace internal
|
||||
|
||||
|
||||
/** If m_matT(i+1,i) is neglegible in floating point arithmetic
|
||||
@@ -268,9 +273,9 @@ std::complex<RealScalar> ei_sqrt(const std::complex<RealScalar> &z)
|
||||
template<typename MatrixType>
|
||||
inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i)
|
||||
{
|
||||
RealScalar d = ei_norm1(m_matT.coeff(i,i)) + ei_norm1(m_matT.coeff(i+1,i+1));
|
||||
RealScalar sd = ei_norm1(m_matT.coeff(i+1,i));
|
||||
if (ei_isMuchSmallerThan(sd, d, NumTraits<RealScalar>::epsilon()))
|
||||
RealScalar d = internal::norm1(m_matT.coeff(i,i)) + internal::norm1(m_matT.coeff(i+1,i+1));
|
||||
RealScalar sd = internal::norm1(m_matT.coeff(i+1,i));
|
||||
if (internal::isMuchSmallerThan(sd, d, NumTraits<RealScalar>::epsilon()))
|
||||
{
|
||||
m_matT.coeffRef(i+1,i) = ComplexScalar(0);
|
||||
return true;
|
||||
@@ -286,7 +291,7 @@ typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::compu
|
||||
if (iter == 10 || iter == 20)
|
||||
{
|
||||
// exceptional shift, taken from http://www.netlib.org/eispack/comqr.f
|
||||
return ei_abs(ei_real(m_matT.coeff(iu,iu-1))) + ei_abs(ei_real(m_matT.coeff(iu-1,iu-2)));
|
||||
return internal::abs(internal::real(m_matT.coeff(iu,iu-1))) + internal::abs(internal::real(m_matT.coeff(iu-1,iu-2)));
|
||||
}
|
||||
|
||||
// compute the shift as one of the eigenvalues of t, the 2x2
|
||||
@@ -297,19 +302,19 @@ typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::compu
|
||||
|
||||
ComplexScalar b = t.coeff(0,1) * t.coeff(1,0);
|
||||
ComplexScalar c = t.coeff(0,0) - t.coeff(1,1);
|
||||
ComplexScalar disc = ei_sqrt(c*c + RealScalar(4)*b);
|
||||
ComplexScalar disc = internal::sqrt(c*c + RealScalar(4)*b);
|
||||
ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b;
|
||||
ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1);
|
||||
ComplexScalar eival1 = (trace + disc) / RealScalar(2);
|
||||
ComplexScalar eival2 = (trace - disc) / RealScalar(2);
|
||||
|
||||
if(ei_norm1(eival1) > ei_norm1(eival2))
|
||||
if(internal::norm1(eival1) > internal::norm1(eival2))
|
||||
eival2 = det / eival1;
|
||||
else
|
||||
eival1 = det / eival2;
|
||||
|
||||
// choose the eigenvalue closest to the bottom entry of the diagonal
|
||||
if(ei_norm1(eival1-t.coeff(1,1)) < ei_norm1(eival2-t.coeff(1,1)))
|
||||
if(internal::norm1(eival1-t.coeff(1,1)) < internal::norm1(eival2-t.coeff(1,1)))
|
||||
return normt * eival1;
|
||||
else
|
||||
return normt * eival2;
|
||||
@@ -320,7 +325,7 @@ template<typename MatrixType>
|
||||
ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
|
||||
{
|
||||
m_matUisUptodate = false;
|
||||
ei_assert(matrix.cols() == matrix.rows());
|
||||
eigen_assert(matrix.cols() == matrix.rows());
|
||||
|
||||
if(matrix.cols() == 1)
|
||||
{
|
||||
@@ -332,14 +337,16 @@ ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& ma
|
||||
return *this;
|
||||
}
|
||||
|
||||
ei_complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix, computeU);
|
||||
internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix, computeU);
|
||||
reduceToTriangularForm(computeU);
|
||||
return *this;
|
||||
}
|
||||
|
||||
namespace internal {
|
||||
|
||||
/* Reduce given matrix to Hessenberg form */
|
||||
template<typename MatrixType, bool IsComplex>
|
||||
struct ei_complex_schur_reduce_to_hessenberg
|
||||
struct complex_schur_reduce_to_hessenberg
|
||||
{
|
||||
// this is the implementation for the case IsComplex = true
|
||||
static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
|
||||
@@ -351,7 +358,7 @@ struct ei_complex_schur_reduce_to_hessenberg
|
||||
};
|
||||
|
||||
template<typename MatrixType>
|
||||
struct ei_complex_schur_reduce_to_hessenberg<MatrixType, false>
|
||||
struct complex_schur_reduce_to_hessenberg<MatrixType, false>
|
||||
{
|
||||
static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
|
||||
{
|
||||
@@ -370,6 +377,8 @@ struct ei_complex_schur_reduce_to_hessenberg<MatrixType, false>
|
||||
}
|
||||
};
|
||||
|
||||
} // end namespace internal
|
||||
|
||||
// Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
|
||||
template<typename MatrixType>
|
||||
void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
|
||||
|
||||
@@ -211,8 +211,8 @@ template<typename _MatrixType> class EigenSolver
|
||||
*/
|
||||
const MatrixType& pseudoEigenvectors() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "EigenSolver is not initialized.");
|
||||
ei_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
|
||||
eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
|
||||
eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
|
||||
return m_eivec;
|
||||
}
|
||||
|
||||
@@ -254,7 +254,7 @@ template<typename _MatrixType> class EigenSolver
|
||||
*/
|
||||
const EigenvalueType& eigenvalues() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "EigenSolver is not initialized.");
|
||||
eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
|
||||
return m_eivalues;
|
||||
}
|
||||
|
||||
@@ -289,7 +289,7 @@ template<typename _MatrixType> class EigenSolver
|
||||
|
||||
ComputationInfo info() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
|
||||
eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
|
||||
return m_realSchur.info();
|
||||
}
|
||||
|
||||
@@ -311,17 +311,17 @@ template<typename _MatrixType> class EigenSolver
|
||||
template<typename MatrixType>
|
||||
MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "EigenSolver is not initialized.");
|
||||
eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
|
||||
Index n = m_eivalues.rows();
|
||||
MatrixType matD = MatrixType::Zero(n,n);
|
||||
for (Index i=0; i<n; ++i)
|
||||
{
|
||||
if (ei_isMuchSmallerThan(ei_imag(m_eivalues.coeff(i)), ei_real(m_eivalues.coeff(i))))
|
||||
matD.coeffRef(i,i) = ei_real(m_eivalues.coeff(i));
|
||||
if (internal::isMuchSmallerThan(internal::imag(m_eivalues.coeff(i)), internal::real(m_eivalues.coeff(i))))
|
||||
matD.coeffRef(i,i) = internal::real(m_eivalues.coeff(i));
|
||||
else
|
||||
{
|
||||
matD.template block<2,2>(i,i) << ei_real(m_eivalues.coeff(i)), ei_imag(m_eivalues.coeff(i)),
|
||||
-ei_imag(m_eivalues.coeff(i)), ei_real(m_eivalues.coeff(i));
|
||||
matD.template block<2,2>(i,i) << internal::real(m_eivalues.coeff(i)), internal::imag(m_eivalues.coeff(i)),
|
||||
-internal::imag(m_eivalues.coeff(i)), internal::real(m_eivalues.coeff(i));
|
||||
++i;
|
||||
}
|
||||
}
|
||||
@@ -331,13 +331,13 @@ MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
|
||||
template<typename MatrixType>
|
||||
typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eigenvectors() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "EigenSolver is not initialized.");
|
||||
ei_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
|
||||
eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
|
||||
eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
|
||||
Index n = m_eivec.cols();
|
||||
EigenvectorsType matV(n,n);
|
||||
for (Index j=0; j<n; ++j)
|
||||
{
|
||||
if (ei_isMuchSmallerThan(ei_imag(m_eivalues.coeff(j)), ei_real(m_eivalues.coeff(j))))
|
||||
if (internal::isMuchSmallerThan(internal::imag(m_eivalues.coeff(j)), internal::real(m_eivalues.coeff(j))))
|
||||
{
|
||||
// we have a real eigen value
|
||||
matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
|
||||
@@ -384,7 +384,7 @@ EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matr
|
||||
else
|
||||
{
|
||||
Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
|
||||
Scalar z = ei_sqrt(ei_abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
|
||||
Scalar z = internal::sqrt(internal::abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
|
||||
m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
|
||||
m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
|
||||
i += 2;
|
||||
@@ -407,7 +407,7 @@ template<typename Scalar>
|
||||
std::complex<Scalar> cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi)
|
||||
{
|
||||
Scalar r,d;
|
||||
if (ei_abs(yr) > ei_abs(yi))
|
||||
if (internal::abs(yr) > internal::abs(yi))
|
||||
{
|
||||
r = yi/yr;
|
||||
d = yr + r*yi;
|
||||
@@ -480,14 +480,14 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
|
||||
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;
|
||||
m_matT.coeffRef(i,n) = t;
|
||||
if (ei_abs(x) > ei_abs(lastw))
|
||||
if (internal::abs(x) > internal::abs(lastw))
|
||||
m_matT.coeffRef(i+1,n) = (-r - w * t) / x;
|
||||
else
|
||||
m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw;
|
||||
}
|
||||
|
||||
// Overflow control
|
||||
Scalar t = ei_abs(m_matT.coeff(i,n));
|
||||
Scalar t = internal::abs(m_matT.coeff(i,n));
|
||||
if ((eps * t) * t > 1)
|
||||
m_matT.col(n).tail(size-i) /= t;
|
||||
}
|
||||
@@ -499,16 +499,16 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
|
||||
Index l = n-1;
|
||||
|
||||
// Last vector component imaginary so matrix is triangular
|
||||
if (ei_abs(m_matT.coeff(n,n-1)) > ei_abs(m_matT.coeff(n-1,n)))
|
||||
if (internal::abs(m_matT.coeff(n,n-1)) > internal::abs(m_matT.coeff(n-1,n)))
|
||||
{
|
||||
m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1);
|
||||
m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1);
|
||||
}
|
||||
else
|
||||
{
|
||||
std::complex<Scalar> cc = cdiv<Scalar>(0.0,-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-p,q);
|
||||
m_matT.coeffRef(n-1,n-1) = ei_real(cc);
|
||||
m_matT.coeffRef(n-1,n) = ei_imag(cc);
|
||||
std::complex<Scalar> cc = cdiv<Scalar>(0.0,-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-p,q);
|
||||
m_matT.coeffRef(n-1,n-1) = internal::real(cc);
|
||||
m_matT.coeffRef(n-1,n) = internal::imag(cc);
|
||||
}
|
||||
m_matT.coeffRef(n,n-1) = 0.0;
|
||||
m_matT.coeffRef(n,n) = 1.0;
|
||||
@@ -530,8 +530,8 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
|
||||
if (m_eivalues.coeff(i).imag() == 0)
|
||||
{
|
||||
std::complex<Scalar> cc = cdiv(-ra,-sa,w,q);
|
||||
m_matT.coeffRef(i,n-1) = ei_real(cc);
|
||||
m_matT.coeffRef(i,n) = ei_imag(cc);
|
||||
m_matT.coeffRef(i,n-1) = internal::real(cc);
|
||||
m_matT.coeffRef(i,n) = internal::imag(cc);
|
||||
}
|
||||
else
|
||||
{
|
||||
@@ -541,12 +541,12 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
|
||||
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(lastw));
|
||||
vr = eps * norm * (internal::abs(w) + internal::abs(q) + internal::abs(x) + internal::abs(y) + internal::abs(lastw));
|
||||
|
||||
std::complex<Scalar> cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi);
|
||||
m_matT.coeffRef(i,n-1) = ei_real(cc);
|
||||
m_matT.coeffRef(i,n) = ei_imag(cc);
|
||||
if (ei_abs(x) > (ei_abs(lastw) + ei_abs(q)))
|
||||
m_matT.coeffRef(i,n-1) = internal::real(cc);
|
||||
m_matT.coeffRef(i,n) = internal::imag(cc);
|
||||
if (internal::abs(x) > (internal::abs(lastw) + internal::abs(q)))
|
||||
{
|
||||
m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x;
|
||||
m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x;
|
||||
@@ -554,13 +554,13 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
|
||||
else
|
||||
{
|
||||
cc = cdiv(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n),lastw,q);
|
||||
m_matT.coeffRef(i+1,n-1) = ei_real(cc);
|
||||
m_matT.coeffRef(i+1,n) = ei_imag(cc);
|
||||
m_matT.coeffRef(i+1,n-1) = internal::real(cc);
|
||||
m_matT.coeffRef(i+1,n) = internal::imag(cc);
|
||||
}
|
||||
}
|
||||
|
||||
// Overflow control
|
||||
Scalar t = std::max(ei_abs(m_matT.coeff(i,n-1)),ei_abs(m_matT.coeff(i,n)));
|
||||
Scalar t = std::max(internal::abs(m_matT.coeff(i,n-1)),internal::abs(m_matT.coeff(i,n)));
|
||||
if ((eps * t) * t > 1)
|
||||
m_matT.block(i, n-1, size-i, 2) /= t;
|
||||
|
||||
|
||||
@@ -182,8 +182,8 @@ template<typename MatrixType>
|
||||
GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>::
|
||||
compute(const MatrixType& matA, const MatrixType& matB, int options)
|
||||
{
|
||||
ei_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
|
||||
ei_assert((options&~(EigVecMask|GenEigMask))==0
|
||||
eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
|
||||
eigen_assert((options&~(EigVecMask|GenEigMask))==0
|
||||
&& (options&EigVecMask)!=EigVecMask
|
||||
&& ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx
|
||||
|| (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx)
|
||||
|
||||
@@ -28,11 +28,13 @@
|
||||
|
||||
template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType;
|
||||
|
||||
namespace internal {
|
||||
template<typename MatrixType>
|
||||
struct ei_traits<HessenbergDecompositionMatrixHReturnType<MatrixType> >
|
||||
struct traits<HessenbergDecompositionMatrixHReturnType<MatrixType> >
|
||||
{
|
||||
typedef MatrixType ReturnType;
|
||||
};
|
||||
}
|
||||
|
||||
/** \eigenvalues_module \ingroup Eigenvalues_Module
|
||||
*
|
||||
@@ -184,7 +186,7 @@ template<typename _MatrixType> class HessenbergDecomposition
|
||||
*/
|
||||
const CoeffVectorType& householderCoefficients() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
|
||||
eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
|
||||
return m_hCoeffs;
|
||||
}
|
||||
|
||||
@@ -219,7 +221,7 @@ template<typename _MatrixType> class HessenbergDecomposition
|
||||
*/
|
||||
const MatrixType& packedMatrix() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
|
||||
eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
|
||||
return m_matrix;
|
||||
}
|
||||
|
||||
@@ -239,7 +241,7 @@ template<typename _MatrixType> class HessenbergDecomposition
|
||||
*/
|
||||
HouseholderSequenceType matrixQ() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
|
||||
eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
|
||||
return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate(), false, m_matrix.rows() - 1, 1);
|
||||
}
|
||||
|
||||
@@ -265,7 +267,7 @@ template<typename _MatrixType> class HessenbergDecomposition
|
||||
*/
|
||||
HessenbergDecompositionMatrixHReturnType<MatrixType> matrixH() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
|
||||
eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
|
||||
return HessenbergDecompositionMatrixHReturnType<MatrixType>(*this);
|
||||
}
|
||||
|
||||
@@ -319,7 +321,7 @@ void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVector
|
||||
|
||||
// A = A H'
|
||||
matA.rightCols(remainingSize)
|
||||
.applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1).conjugate(), ei_conj(h), &temp.coeffRef(0));
|
||||
.applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1).conjugate(), internal::conj(h), &temp.coeffRef(0));
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@@ -26,10 +26,10 @@
|
||||
#ifndef EIGEN_MATRIXBASEEIGENVALUES_H
|
||||
#define EIGEN_MATRIXBASEEIGENVALUES_H
|
||||
|
||||
|
||||
namespace internal {
|
||||
|
||||
template<typename Derived, bool IsComplex>
|
||||
struct ei_eigenvalues_selector
|
||||
struct eigenvalues_selector
|
||||
{
|
||||
// this is the implementation for the case IsComplex = true
|
||||
static inline typename MatrixBase<Derived>::EigenvaluesReturnType const
|
||||
@@ -42,7 +42,7 @@ struct ei_eigenvalues_selector
|
||||
};
|
||||
|
||||
template<typename Derived>
|
||||
struct ei_eigenvalues_selector<Derived, false>
|
||||
struct eigenvalues_selector<Derived, false>
|
||||
{
|
||||
static inline typename MatrixBase<Derived>::EigenvaluesReturnType const
|
||||
run(const MatrixBase<Derived>& m)
|
||||
@@ -53,6 +53,8 @@ struct ei_eigenvalues_selector<Derived, false>
|
||||
}
|
||||
};
|
||||
|
||||
} // end namespace internal
|
||||
|
||||
/** \brief Computes the eigenvalues of a matrix
|
||||
* \returns Column vector containing the eigenvalues.
|
||||
*
|
||||
@@ -77,8 +79,8 @@ template<typename Derived>
|
||||
inline typename MatrixBase<Derived>::EigenvaluesReturnType
|
||||
MatrixBase<Derived>::eigenvalues() const
|
||||
{
|
||||
typedef typename ei_traits<Derived>::Scalar Scalar;
|
||||
return ei_eigenvalues_selector<Derived, NumTraits<Scalar>::IsComplex>::run(derived());
|
||||
typedef typename internal::traits<Derived>::Scalar Scalar;
|
||||
return internal::eigenvalues_selector<Derived, NumTraits<Scalar>::IsComplex>::run(derived());
|
||||
}
|
||||
|
||||
/** \brief Computes the eigenvalues of a matrix
|
||||
@@ -135,7 +137,7 @@ MatrixBase<Derived>::operatorNorm() const
|
||||
typename Derived::PlainObject m_eval(derived());
|
||||
// FIXME if it is really guaranteed that the eigenvalues are already sorted,
|
||||
// then we don't need to compute a maxCoeff() here, comparing the 1st and last ones is enough.
|
||||
return ei_sqrt((m_eval*m_eval.adjoint())
|
||||
return internal::sqrt((m_eval*m_eval.adjoint())
|
||||
.eval()
|
||||
.template selfadjointView<Lower>()
|
||||
.eigenvalues()
|
||||
|
||||
@@ -137,8 +137,8 @@ template<typename _MatrixType> class RealSchur
|
||||
*/
|
||||
const MatrixType& matrixU() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "RealSchur is not initialized.");
|
||||
ei_assert(m_matUisUptodate && "The matrix U has not been computed during the RealSchur decomposition.");
|
||||
eigen_assert(m_isInitialized && "RealSchur is not initialized.");
|
||||
eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the RealSchur decomposition.");
|
||||
return m_matU;
|
||||
}
|
||||
|
||||
@@ -154,7 +154,7 @@ template<typename _MatrixType> class RealSchur
|
||||
*/
|
||||
const MatrixType& matrixT() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "RealSchur is not initialized.");
|
||||
eigen_assert(m_isInitialized && "RealSchur is not initialized.");
|
||||
return m_matT;
|
||||
}
|
||||
|
||||
@@ -183,7 +183,7 @@ template<typename _MatrixType> class RealSchur
|
||||
*/
|
||||
ComputationInfo info() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "RealSchur is not initialized.");
|
||||
eigen_assert(m_isInitialized && "RealSchur is not initialized.");
|
||||
return m_info;
|
||||
}
|
||||
|
||||
@@ -300,10 +300,10 @@ inline typename MatrixType::Index RealSchur<MatrixType>::findSmallSubdiagEntry(I
|
||||
Index res = iu;
|
||||
while (res > 0)
|
||||
{
|
||||
Scalar s = ei_abs(m_matT.coeff(res-1,res-1)) + ei_abs(m_matT.coeff(res,res));
|
||||
Scalar s = internal::abs(m_matT.coeff(res-1,res-1)) + internal::abs(m_matT.coeff(res,res));
|
||||
if (s == 0.0)
|
||||
s = norm;
|
||||
if (ei_abs(m_matT.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
|
||||
if (internal::abs(m_matT.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
|
||||
break;
|
||||
res--;
|
||||
}
|
||||
@@ -325,7 +325,7 @@ inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, Scal
|
||||
|
||||
if (q >= 0) // Two real eigenvalues
|
||||
{
|
||||
Scalar z = ei_sqrt(ei_abs(q));
|
||||
Scalar z = internal::sqrt(internal::abs(q));
|
||||
JacobiRotation<Scalar> rot;
|
||||
if (p >= 0)
|
||||
rot.makeGivens(p + z, m_matT.coeff(iu, iu-1));
|
||||
@@ -357,7 +357,7 @@ inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& ex
|
||||
exshift += shiftInfo.coeff(0);
|
||||
for (Index i = 0; i <= iu; ++i)
|
||||
m_matT.coeffRef(i,i) -= shiftInfo.coeff(0);
|
||||
Scalar s = ei_abs(m_matT.coeff(iu,iu-1)) + ei_abs(m_matT.coeff(iu-1,iu-2));
|
||||
Scalar s = internal::abs(m_matT.coeff(iu,iu-1)) + internal::abs(m_matT.coeff(iu-1,iu-2));
|
||||
shiftInfo.coeffRef(0) = Scalar(0.75) * s;
|
||||
shiftInfo.coeffRef(1) = Scalar(0.75) * s;
|
||||
shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s;
|
||||
@@ -370,7 +370,7 @@ inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& ex
|
||||
s = s * s + shiftInfo.coeff(2);
|
||||
if (s > 0)
|
||||
{
|
||||
s = ei_sqrt(s);
|
||||
s = internal::sqrt(s);
|
||||
if (shiftInfo.coeff(1) < shiftInfo.coeff(0))
|
||||
s = -s;
|
||||
s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
|
||||
@@ -400,9 +400,9 @@ inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const V
|
||||
if (im == il) {
|
||||
break;
|
||||
}
|
||||
const Scalar lhs = m_matT.coeff(im,im-1) * (ei_abs(v.coeff(1)) + ei_abs(v.coeff(2)));
|
||||
const Scalar rhs = v.coeff(0) * (ei_abs(m_matT.coeff(im-1,im-1)) + ei_abs(Tmm) + ei_abs(m_matT.coeff(im+1,im+1)));
|
||||
if (ei_abs(lhs) < NumTraits<Scalar>::epsilon() * rhs)
|
||||
const Scalar lhs = m_matT.coeff(im,im-1) * (internal::abs(v.coeff(1)) + internal::abs(v.coeff(2)));
|
||||
const Scalar rhs = v.coeff(0) * (internal::abs(m_matT.coeff(im-1,im-1)) + internal::abs(Tmm) + internal::abs(m_matT.coeff(im+1,im+1)));
|
||||
if (internal::abs(lhs) < NumTraits<Scalar>::epsilon() * rhs)
|
||||
{
|
||||
break;
|
||||
}
|
||||
|
||||
@@ -101,7 +101,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
|
||||
* This is a column vector with entries of type #RealScalar.
|
||||
* The length of the vector is the size of \p _MatrixType.
|
||||
*/
|
||||
typedef typename ei_plain_col_type<MatrixType, RealScalar>::type RealVectorType;
|
||||
typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType;
|
||||
typedef Tridiagonalization<MatrixType> TridiagonalizationType;
|
||||
|
||||
/** \brief Default constructor for fixed-size matrices.
|
||||
@@ -219,8 +219,8 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
|
||||
*/
|
||||
const MatrixType& eigenvectors() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
|
||||
ei_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
|
||||
eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
|
||||
eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
|
||||
return m_eivec;
|
||||
}
|
||||
|
||||
@@ -240,7 +240,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
|
||||
*/
|
||||
const RealVectorType& eigenvalues() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
|
||||
eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
|
||||
return m_eivalues;
|
||||
}
|
||||
|
||||
@@ -264,8 +264,8 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
|
||||
*/
|
||||
MatrixType operatorSqrt() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
|
||||
ei_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
|
||||
eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
|
||||
eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
|
||||
return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
|
||||
}
|
||||
|
||||
@@ -289,8 +289,8 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
|
||||
*/
|
||||
MatrixType operatorInverseSqrt() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
|
||||
ei_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
|
||||
eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
|
||||
eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
|
||||
return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
|
||||
}
|
||||
|
||||
@@ -300,7 +300,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
|
||||
*/
|
||||
ComputationInfo info() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
|
||||
eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
|
||||
return m_info;
|
||||
}
|
||||
|
||||
@@ -335,15 +335,17 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
|
||||
* Implemented from Golub's "Matrix Computations", algorithm 8.3.2:
|
||||
* "implicit symmetric QR step with Wilkinson shift"
|
||||
*/
|
||||
namespace internal {
|
||||
template<typename RealScalar, typename Scalar, typename Index>
|
||||
static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n);
|
||||
static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n);
|
||||
}
|
||||
|
||||
template<typename MatrixType>
|
||||
SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
|
||||
::compute(const MatrixType& matrix, int options)
|
||||
{
|
||||
ei_assert(matrix.cols() == matrix.rows());
|
||||
ei_assert((options&~(EigVecMask|GenEigMask))==0
|
||||
eigen_assert(matrix.cols() == matrix.rows());
|
||||
eigen_assert((options&~(EigVecMask|GenEigMask))==0
|
||||
&& (options&EigVecMask)!=EigVecMask
|
||||
&& "invalid option parameter");
|
||||
bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
|
||||
@@ -352,7 +354,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
|
||||
|
||||
if(n==1)
|
||||
{
|
||||
m_eivalues.coeffRef(0,0) = ei_real(matrix.coeff(0,0));
|
||||
m_eivalues.coeffRef(0,0) = internal::real(matrix.coeff(0,0));
|
||||
if(computeEigenvectors)
|
||||
m_eivec.setOnes();
|
||||
m_info = Success;
|
||||
@@ -367,7 +369,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
|
||||
|
||||
mat = matrix;
|
||||
m_subdiag.resize(n-1);
|
||||
ei_tridiagonalization_inplace(mat, diag, m_subdiag, computeEigenvectors);
|
||||
internal::tridiagonalization_inplace(mat, diag, m_subdiag, computeEigenvectors);
|
||||
|
||||
Index end = n-1;
|
||||
Index start = 0;
|
||||
@@ -376,7 +378,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
|
||||
while (end>0)
|
||||
{
|
||||
for (Index i = start; i<end; ++i)
|
||||
if (ei_isMuchSmallerThan(ei_abs(m_subdiag[i]),(ei_abs(diag[i])+ei_abs(diag[i+1]))))
|
||||
if (internal::isMuchSmallerThan(internal::abs(m_subdiag[i]),(internal::abs(diag[i])+internal::abs(diag[i+1]))))
|
||||
m_subdiag[i] = 0;
|
||||
|
||||
// find the largest unreduced block
|
||||
@@ -396,7 +398,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
|
||||
while (start>0 && m_subdiag[start-1]!=0)
|
||||
start--;
|
||||
|
||||
ei_tridiagonal_qr_step(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n);
|
||||
internal::tridiagonal_qr_step(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n);
|
||||
}
|
||||
|
||||
if (iter <= m_maxIterations)
|
||||
@@ -427,12 +429,13 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
|
||||
return *this;
|
||||
}
|
||||
|
||||
namespace internal {
|
||||
template<typename RealScalar, typename Scalar, typename Index>
|
||||
static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n)
|
||||
static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n)
|
||||
{
|
||||
RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
|
||||
RealScalar e2 = ei_abs2(subdiag[end-1]);
|
||||
RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * ei_sqrt(td*td + e2));
|
||||
RealScalar e2 = abs2(subdiag[end-1]);
|
||||
RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2));
|
||||
RealScalar x = diag[start] - mu;
|
||||
RealScalar z = subdiag[start];
|
||||
|
||||
@@ -468,5 +471,6 @@ static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index
|
||||
}
|
||||
}
|
||||
}
|
||||
} // end namespace internal
|
||||
|
||||
#endif // EIGEN_SELFADJOINTEIGENSOLVER_H
|
||||
|
||||
@@ -26,6 +26,11 @@
|
||||
#ifndef EIGEN_TRIDIAGONALIZATION_H
|
||||
#define EIGEN_TRIDIAGONALIZATION_H
|
||||
|
||||
namespace internal {
|
||||
template<typename MatrixType, typename CoeffVectorType>
|
||||
void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs);
|
||||
}
|
||||
|
||||
/** \eigenvalues_module \ingroup Eigenvalues_Module
|
||||
*
|
||||
*
|
||||
@@ -78,15 +83,15 @@ template<typename _MatrixType> class Tridiagonalization
|
||||
};
|
||||
|
||||
typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
|
||||
typedef typename ei_plain_col_type<MatrixType, RealScalar>::type DiagonalType;
|
||||
typedef typename internal::plain_col_type<MatrixType, RealScalar>::type DiagonalType;
|
||||
typedef Matrix<RealScalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> SubDiagonalType;
|
||||
|
||||
typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex,
|
||||
typedef typename internal::meta_if<NumTraits<Scalar>::IsComplex,
|
||||
typename Diagonal<MatrixType,0>::RealReturnType,
|
||||
Diagonal<MatrixType,0>
|
||||
>::ret DiagonalReturnType;
|
||||
|
||||
typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex,
|
||||
typedef typename internal::meta_if<NumTraits<Scalar>::IsComplex,
|
||||
typename Diagonal<
|
||||
Block<MatrixType,SizeMinusOne,SizeMinusOne>,0 >::RealReturnType,
|
||||
Diagonal<
|
||||
@@ -129,7 +134,7 @@ template<typename _MatrixType> class Tridiagonalization
|
||||
m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1),
|
||||
m_isInitialized(false)
|
||||
{
|
||||
ei_tridiagonalization_inplace(m_matrix, m_hCoeffs);
|
||||
internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
|
||||
m_isInitialized = true;
|
||||
}
|
||||
|
||||
@@ -154,7 +159,7 @@ template<typename _MatrixType> class Tridiagonalization
|
||||
{
|
||||
m_matrix = matrix;
|
||||
m_hCoeffs.resize(matrix.rows()-1, 1);
|
||||
ei_tridiagonalization_inplace(m_matrix, m_hCoeffs);
|
||||
internal::tridiagonalization_inplace(m_matrix, m_hCoeffs);
|
||||
m_isInitialized = true;
|
||||
return *this;
|
||||
}
|
||||
@@ -177,7 +182,7 @@ template<typename _MatrixType> class Tridiagonalization
|
||||
*/
|
||||
inline CoeffVectorType householderCoefficients() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "Tridiagonalization is not initialized.");
|
||||
eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
|
||||
return m_hCoeffs;
|
||||
}
|
||||
|
||||
@@ -214,7 +219,7 @@ template<typename _MatrixType> class Tridiagonalization
|
||||
*/
|
||||
inline const MatrixType& packedMatrix() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "Tridiagonalization is not initialized.");
|
||||
eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
|
||||
return m_matrix;
|
||||
}
|
||||
|
||||
@@ -235,7 +240,7 @@ template<typename _MatrixType> class Tridiagonalization
|
||||
*/
|
||||
HouseholderSequenceType matrixQ() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "Tridiagonalization is not initialized.");
|
||||
eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
|
||||
return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate(), false, m_matrix.rows() - 1, 1);
|
||||
}
|
||||
|
||||
@@ -296,7 +301,7 @@ template<typename MatrixType>
|
||||
const typename Tridiagonalization<MatrixType>::DiagonalReturnType
|
||||
Tridiagonalization<MatrixType>::diagonal() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "Tridiagonalization is not initialized.");
|
||||
eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
|
||||
return m_matrix.diagonal();
|
||||
}
|
||||
|
||||
@@ -304,7 +309,7 @@ template<typename MatrixType>
|
||||
const typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
|
||||
Tridiagonalization<MatrixType>::subDiagonal() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "Tridiagonalization is not initialized.");
|
||||
eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
|
||||
Index n = m_matrix.rows();
|
||||
return Block<MatrixType,SizeMinusOne,SizeMinusOne>(m_matrix, 1, 0, n-1,n-1).diagonal();
|
||||
}
|
||||
@@ -315,7 +320,7 @@ Tridiagonalization<MatrixType>::matrixT() const
|
||||
{
|
||||
// FIXME should this function (and other similar ones) rather take a matrix as argument
|
||||
// and fill it ? (to avoid temporaries)
|
||||
ei_assert(m_isInitialized && "Tridiagonalization is not initialized.");
|
||||
eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
|
||||
Index n = m_matrix.rows();
|
||||
MatrixType matT = m_matrix;
|
||||
matT.topRightCorner(n-1, n-1).diagonal() = subDiagonal().template cast<Scalar>().conjugate();
|
||||
@@ -327,6 +332,8 @@ Tridiagonalization<MatrixType>::matrixT() const
|
||||
return matT;
|
||||
}
|
||||
|
||||
namespace internal {
|
||||
|
||||
/** \internal
|
||||
* Performs a tridiagonal decomposition of the selfadjoint matrix \a matA in-place.
|
||||
*
|
||||
@@ -351,10 +358,10 @@ Tridiagonalization<MatrixType>::matrixT() const
|
||||
* \sa Tridiagonalization::packedMatrix()
|
||||
*/
|
||||
template<typename MatrixType, typename CoeffVectorType>
|
||||
void ei_tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs)
|
||||
void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs)
|
||||
{
|
||||
ei_assert(matA.rows()==matA.cols());
|
||||
ei_assert(matA.rows()==hCoeffs.size()+1);
|
||||
eigen_assert(matA.rows()==matA.cols());
|
||||
eigen_assert(matA.rows()==hCoeffs.size()+1);
|
||||
typedef typename MatrixType::Index Index;
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::RealScalar RealScalar;
|
||||
@@ -371,9 +378,9 @@ void ei_tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs)
|
||||
matA.col(i).coeffRef(i+1) = 1;
|
||||
|
||||
hCoeffs.tail(n-i-1).noalias() = (matA.bottomRightCorner(remainingSize,remainingSize).template selfadjointView<Lower>()
|
||||
* (ei_conj(h) * matA.col(i).tail(remainingSize)));
|
||||
* (conj(h) * matA.col(i).tail(remainingSize)));
|
||||
|
||||
hCoeffs.tail(n-i-1) += (ei_conj(h)*Scalar(-0.5)*(hCoeffs.tail(remainingSize).dot(matA.col(i).tail(remainingSize)))) * matA.col(i).tail(n-i-1);
|
||||
hCoeffs.tail(n-i-1) += (conj(h)*Scalar(-0.5)*(hCoeffs.tail(remainingSize).dot(matA.col(i).tail(remainingSize)))) * matA.col(i).tail(n-i-1);
|
||||
|
||||
matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>()
|
||||
.rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), -1);
|
||||
@@ -387,7 +394,7 @@ void ei_tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs)
|
||||
template<typename MatrixType,
|
||||
int Size=MatrixType::ColsAtCompileTime,
|
||||
bool IsComplex=NumTraits<typename MatrixType::Scalar>::IsComplex>
|
||||
struct ei_tridiagonalization_inplace_selector;
|
||||
struct tridiagonalization_inplace_selector;
|
||||
|
||||
/** \brief Performs a full tridiagonalization in place
|
||||
*
|
||||
@@ -430,19 +437,19 @@ struct ei_tridiagonalization_inplace_selector;
|
||||
* \sa class Tridiagonalization
|
||||
*/
|
||||
template<typename MatrixType, typename DiagonalType, typename SubDiagonalType>
|
||||
void ei_tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
|
||||
void tridiagonalization_inplace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
|
||||
{
|
||||
typedef typename MatrixType::Index Index;
|
||||
//Index n = mat.rows();
|
||||
ei_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1);
|
||||
ei_tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, extractQ);
|
||||
eigen_assert(mat.cols()==mat.rows() && diag.size()==mat.rows() && subdiag.size()==mat.rows()-1);
|
||||
tridiagonalization_inplace_selector<MatrixType>::run(mat, diag, subdiag, extractQ);
|
||||
}
|
||||
|
||||
/** \internal
|
||||
* General full tridiagonalization
|
||||
*/
|
||||
template<typename MatrixType, int Size, bool IsComplex>
|
||||
struct ei_tridiagonalization_inplace_selector
|
||||
struct tridiagonalization_inplace_selector
|
||||
{
|
||||
typedef typename Tridiagonalization<MatrixType>::CoeffVectorType CoeffVectorType;
|
||||
typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType;
|
||||
@@ -451,7 +458,7 @@ struct ei_tridiagonalization_inplace_selector
|
||||
static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
|
||||
{
|
||||
CoeffVectorType hCoeffs(mat.cols()-1);
|
||||
ei_tridiagonalization_inplace(mat,hCoeffs);
|
||||
tridiagonalization_inplace(mat,hCoeffs);
|
||||
diag = mat.diagonal().real();
|
||||
subdiag = mat.template diagonal<-1>().real();
|
||||
if(extractQ)
|
||||
@@ -464,7 +471,7 @@ struct ei_tridiagonalization_inplace_selector
|
||||
* Especially useful for plane fitting.
|
||||
*/
|
||||
template<typename MatrixType>
|
||||
struct ei_tridiagonalization_inplace_selector<MatrixType,3,false>
|
||||
struct tridiagonalization_inplace_selector<MatrixType,3,false>
|
||||
{
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::RealScalar RealScalar;
|
||||
@@ -473,7 +480,7 @@ struct ei_tridiagonalization_inplace_selector<MatrixType,3,false>
|
||||
static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
|
||||
{
|
||||
diag[0] = mat(0,0);
|
||||
RealScalar v1norm2 = ei_abs2(mat(2,0));
|
||||
RealScalar v1norm2 = abs2(mat(2,0));
|
||||
if(v1norm2 == RealScalar(0))
|
||||
{
|
||||
diag[1] = mat(1,1);
|
||||
@@ -485,7 +492,7 @@ struct ei_tridiagonalization_inplace_selector<MatrixType,3,false>
|
||||
}
|
||||
else
|
||||
{
|
||||
RealScalar beta = ei_sqrt(ei_abs2(mat(1,0)) + v1norm2);
|
||||
RealScalar beta = sqrt(abs2(mat(1,0)) + v1norm2);
|
||||
RealScalar invBeta = RealScalar(1)/beta;
|
||||
Scalar m01 = mat(1,0) * invBeta;
|
||||
Scalar m02 = mat(2,0) * invBeta;
|
||||
@@ -508,16 +515,19 @@ struct ei_tridiagonalization_inplace_selector<MatrixType,3,false>
|
||||
* Trivial specialization for 1x1 matrices
|
||||
*/
|
||||
template<typename MatrixType, bool IsComplex>
|
||||
struct ei_tridiagonalization_inplace_selector<MatrixType,1,IsComplex>
|
||||
struct tridiagonalization_inplace_selector<MatrixType,1,IsComplex>
|
||||
{
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
|
||||
template<typename DiagonalType, typename SubDiagonalType>
|
||||
static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType&, bool extractQ)
|
||||
{
|
||||
diag(0,0) = ei_real(mat(0,0));
|
||||
diag(0,0) = real(mat(0,0));
|
||||
if(extractQ)
|
||||
mat(0,0) = Scalar(1);
|
||||
}
|
||||
};
|
||||
|
||||
} // end namespace internal
|
||||
|
||||
#endif // EIGEN_TRIDIAGONALIZATION_H
|
||||
|
||||
Reference in New Issue
Block a user