merge and add start/end to Eigen2Support

This commit is contained in:
Gael Guennebaud
2010-01-05 13:07:32 +01:00
78 changed files with 1072 additions and 592 deletions

View File

@@ -133,7 +133,7 @@ void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
for (int i=0; i<n; i++)
{
int k;
m_eivalues.cwiseAbs().end(n-i).minCoeff(&k);
m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k);
if (k != 0)
{
k += i;

View File

@@ -620,7 +620,7 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
// Overflow control
t = ei_abs(matH.coeff(i,n));
if ((eps * t) * t > 1)
matH.col(n).end(nn-i) /= t;
matH.col(n).tail(nn-i) /= t;
}
}
}
@@ -708,7 +708,7 @@ void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
// in this algo low==0 and high==nn-1 !!
if (i < low || i > high)
{
m_eivec.row(i).end(nn-i) = matH.row(i).end(nn-i);
m_eivec.row(i).tail(nn-i) = matH.row(i).tail(nn-i);
}
}

View File

@@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
@@ -55,25 +55,23 @@ template<typename _MatrixType> class HessenbergDecomposition
};
typedef Matrix<Scalar, SizeMinusOne, 1> CoeffVectorType;
typedef Matrix<RealScalar, Size, 1> DiagonalType;
typedef Matrix<RealScalar, SizeMinusOne, 1> SubDiagonalType;
typedef typename Diagonal<MatrixType,0>::RealReturnType DiagonalReturnType;
typedef typename Diagonal<
Block<MatrixType,SizeMinusOne,SizeMinusOne>,0 >::RealReturnType SubDiagonalReturnType;
/** This constructor initializes a HessenbergDecomposition object for
* further use with HessenbergDecomposition::compute()
*/
HessenbergDecomposition(int size = Size==Dynamic ? 2 : Size)
: m_matrix(size,size), m_hCoeffs(size-1)
{}
: m_matrix(size,size)
{
if(size>1)
m_hCoeffs.resize(size-1);
}
HessenbergDecomposition(const MatrixType& matrix)
: m_matrix(matrix),
m_hCoeffs(matrix.cols()-1)
: m_matrix(matrix)
{
if(matrix.rows()<2)
return;
m_hCoeffs.resize(matrix.rows()-1,1);
_compute(m_matrix, m_hCoeffs);
}
@@ -84,6 +82,8 @@ template<typename _MatrixType> class HessenbergDecomposition
void compute(const MatrixType& matrix)
{
m_matrix = matrix;
if(matrix.rows()<2)
return;
m_hCoeffs.resize(matrix.rows()-1,1);
_compute(m_matrix, m_hCoeffs);
}
@@ -150,7 +150,7 @@ void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVector
int remainingSize = n-i-1;
RealScalar beta;
Scalar h;
matA.col(i).end(remainingSize).makeHouseholderInPlace(h, beta);
matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
matA.col(i).coeffRef(i+1) = beta;
hCoeffs.coeffRef(i) = h;
@@ -159,11 +159,11 @@ void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVector
// A = H A
matA.corner(BottomRight, remainingSize, remainingSize)
.applyHouseholderOnTheLeft(matA.col(i).end(remainingSize-1), h, &temp.coeffRef(0));
.applyHouseholderOnTheLeft(matA.col(i).tail(remainingSize-1), h, &temp.coeffRef(0));
// A = A H'
matA.corner(BottomRight, n, remainingSize)
.applyHouseholderOnTheRight(matA.col(i).end(remainingSize-1).conjugate(), ei_conj(h), &temp.coeffRef(0));
.applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1).conjugate(), ei_conj(h), &temp.coeffRef(0));
}
}
@@ -178,7 +178,7 @@ HessenbergDecomposition<MatrixType>::matrixQ() const
for (int i = n-2; i>=0; i--)
{
matQ.corner(BottomRight,n-i-1,n-i-1)
.applyHouseholderOnTheLeft(m_matrix.col(i).end(n-i-2), ei_conj(m_hCoeffs.coeff(i)), &temp.coeffRef(0,0));
.applyHouseholderOnTheLeft(m_matrix.col(i).tail(n-i-2), ei_conj(m_hCoeffs.coeff(i)), &temp.coeffRef(0,0));
}
return matQ;
}

View File

@@ -197,25 +197,24 @@ void Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType&
{
assert(matA.rows()==matA.cols());
int n = matA.rows();
Matrix<Scalar,1,Dynamic> aux(n);
for (int i = 0; i<n-1; ++i)
{
int remainingSize = n-i-1;
RealScalar beta;
Scalar h;
matA.col(i).end(remainingSize).makeHouseholderInPlace(h, beta);
matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
// Apply similarity transformation to remaining columns,
// i.e., A = H A H' where H = I - h v v' and v = matA.col(i).end(n-i-1)
// i.e., A = H A H' where H = I - h v v' and v = matA.col(i).tail(n-i-1)
matA.col(i).coeffRef(i+1) = 1;
hCoeffs.end(n-i-1) = (matA.corner(BottomRight,remainingSize,remainingSize).template selfadjointView<LowerTriangular>()
* (ei_conj(h) * matA.col(i).end(remainingSize)));
hCoeffs.tail(n-i-1) = (matA.corner(BottomRight,remainingSize,remainingSize).template selfadjointView<LowerTriangular>()
* (ei_conj(h) * matA.col(i).tail(remainingSize)));
hCoeffs.end(n-i-1) += (ei_conj(h)*Scalar(-0.5)*(hCoeffs.end(remainingSize).dot(matA.col(i).end(remainingSize)))) * matA.col(i).end(n-i-1);
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);
matA.corner(BottomRight, remainingSize, remainingSize).template selfadjointView<LowerTriangular>()
.rankUpdate(matA.col(i).end(remainingSize), hCoeffs.end(remainingSize), -1);
.rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), -1);
matA.col(i).coeffRef(i+1) = beta;
hCoeffs.coeffRef(i) = h;
@@ -243,7 +242,7 @@ void Tridiagonalization<MatrixType>::matrixQInPlace(MatrixBase<QDerived>* q) con
for (int i = n-2; i>=0; i--)
{
matQ.corner(BottomRight,n-i-1,n-i-1)
.applyHouseholderOnTheLeft(m_matrix.col(i).end(n-i-2), ei_conj(m_hCoeffs.coeff(i)), &aux.coeffRef(0,0));
.applyHouseholderOnTheLeft(m_matrix.col(i).tail(n-i-2), ei_conj(m_hCoeffs.coeff(i)), &aux.coeffRef(0,0));
}
}