mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-18 14:34:17 +08:00
Code Issues Packages Projects Releases Wiki Activity

* added ei_sqrt for complex

Browse Source 
* updated Cholesky to support complex
* correct result_type for abs and abs2 functors
This commit is contained in:
Gael Guennebaud 2008-04-27 14:05:40 +00:00
parent 4ffffa670e
commit 64bacf1c3f
7 changed files with 98 additions and 64 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

44
Eigen/src/Cholesky/Cholesky.h
View File

@ -54,11 +54,6 @@ template<typename MatrixType> class Cholesky
compute(matrix);
}
Triangular<Upper, Temporary<Transpose<MatrixType> > > matrixU(void) const
{
return m_matrix.transpose().temporary().upper();
}
Triangular<Lower, MatrixType> matrixL(void) const
{
return m_matrix.lower();
@ -88,24 +83,9 @@ void Cholesky<MatrixType>::compute(const MatrixType& matrix)
{
assert(matrix.rows()==matrix.cols());
const int size = matrix.rows();
m_matrix = matrix;
m_matrix = matrix.conjugate();
#if 1
// this version looks faster for large matrices
m_isPositiveDefinite = m_matrix(0,0) > Scalar(0);
m_matrix(0,0) = ei_sqrt(m_matrix(0,0));
m_matrix.col(0).end(size-1) = m_matrix.row(0).end(size-1) / m_matrix(0,0);
for (int j = 1; j < size; ++j)
{
Scalar tmp = m_matrix(j,j) - m_matrix.row(j).start(j).norm2();
m_isPositiveDefinite = m_isPositiveDefinite && tmp > Scalar(0);
m_matrix(j,j) = ei_sqrt(tmp<Scalar(0) ? Scalar(0) : tmp);
tmp = Scalar(1) / m_matrix(j,j);
for (int i = j+1; i < size; ++i)
m_matrix(i,j) = tmp * (m_matrix(j,i) -
(m_matrix.row(i).start(j) * m_matrix.row(j).start(j).transpose())(0,0) );
}
#else
#if 0
m_isPositiveDefinite = true;
for (int i = 0; i < size; ++i)
{
@ -118,6 +98,23 @@ void Cholesky<MatrixType>::compute(const MatrixType& matrix)
m_matrix.col(j).end(size-j) -= m_matrix(j,i) * m_matrix.col(i).end(size-j);
}
}
#else
// this version looks faster for large matrices
// m_isPositiveDefinite = m_matrix(0,0) > Scalar(0);
m_matrix(0,0) = ei_sqrt(m_matrix(0,0));
m_matrix.col(0).end(size-1) = m_matrix.row(0).end(size-1) / m_matrix(0,0);
for (int j = 1; j < size; ++j)
{
// Scalar tmp = m_matrix(j,j) - m_matrix.row(j).start(j).norm2();
Scalar tmp = m_matrix(j,j) - (m_matrix.row(j).start(j) * m_matrix.row(j).start(j).adjoint())(0,0);
// m_isPositiveDefinite = m_isPositiveDefinite && tmp > Scalar(0);
// m_matrix(j,j) = ei_sqrt(tmp<Scalar(0) ? Scalar(0) : tmp);
m_matrix(j,j) = ei_sqrt(tmp);
tmp = 1. / m_matrix(j,j);
for (int i = j+1; i < size; ++i)
m_matrix(i,j) = tmp * (m_matrix(j,i) -
(m_matrix.row(i).start(j) * m_matrix.row(j).start(j).adjoint())(0,0) );
}
#endif
}
@ -132,8 +129,7 @@ typename DerivedVec::Eval Cholesky<MatrixType>::solve(MatrixBase<DerivedVec> &ve
// FIXME .inverseProduct creates a temporary that is not nice since it is called twice
// add a .inverseProductInPlace ??
return m_matrix.transpose().upper()
.inverseProduct(m_matrix.lower().inverseProduct(vecB));
return m_matrix.adjoint().upper().inverseProduct(m_matrix.lower().inverseProduct(vecB));
}

21
Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h
View File

@ -55,12 +55,7 @@ template<typename MatrixType> class CholeskyWithoutSquareRoot
compute(matrix);
}
Triangular<Upper|UnitDiagBit, Temporary<Transpose<MatrixType> > > matrixU(void) const
{
return m_matrix.transpose().temporary().upperWithUnitDiag();
}
Triangular<Upper|UnitDiagBit, MatrixType > matrixL(void) const
Triangular<Lower|UnitDiagBit, MatrixType > matrixL(void) const
{
return m_matrix.lowerWithUnitDiag();
}
@ -78,7 +73,7 @@ template<typename MatrixType> class CholeskyWithoutSquareRoot
template<typename DerivedVec>
typename DerivedVec::Eval solve(MatrixBase<DerivedVec> &vecB);
/** Compute the Cholesky decomposition A = U'DU = LDL' of \a matrix
/** Compute / recompute the Cholesky decomposition A = U'DU = LDL' of \a matrix
*/
void compute(const MatrixType& matrix);
@ -97,7 +92,7 @@ void CholeskyWithoutSquareRoot<MatrixType>::compute(const MatrixType& matrix)
{
assert(matrix.rows()==matrix.cols());
const int size = matrix.rows();
m_matrix = matrix;
m_matrix = matrix.conjugate();
#if 0
for (int i = 0; i < size; ++i)
{
@ -113,12 +108,12 @@ void CholeskyWithoutSquareRoot<MatrixType>::compute(const MatrixType& matrix)
m_matrix.col(0).end(size-1) = m_matrix.row(0).end(size-1) / m_matrix(0,0);
for (int j = 1; j < size; ++j)
{
Scalar tmp = m_matrix(j,j) - (m_matrix.row(j).start(j) * m_matrix.col(j).start(j))(0,0);
Scalar tmp = m_matrix(j,j) - (m_matrix.row(j).start(j) * m_matrix.col(j).start(j).conjugate())(0,0);
m_matrix(j,j) = tmp;
tmp = Scalar(1) / tmp;
for (int i = j+1; i < size; ++i)
{
m_matrix(j,i) = (m_matrix(j,i) - (m_matrix.row(i).start(j) * m_matrix.col(j).start(j))(0,0) );
m_matrix(j,i) = (m_matrix(j,i) - (m_matrix.row(i).start(j) * m_matrix.col(j).start(j).conjugate())(0,0) );
m_matrix(i,j) = tmp * m_matrix(j,i);
}
}
@ -136,12 +131,16 @@ typename DerivedVec::Eval CholeskyWithoutSquareRoot<MatrixType>::solve(MatrixBas
// FIXME .inverseProduct creates a temporary that is not nice since it is called twice
// maybe add a .inverseProductInPlace() ??
return m_matrix.transpose().upperWithUnitDiag()
return m_matrix.adjoint().upperWithUnitDiag()
.inverseProduct(
(m_matrix.lowerWithUnitDiag()
.inverseProduct(vecB))
.cwiseQuotient(m_matrix.diagonal())
);
// return m_matrix.adjoint().upperWithUnitDiag()
// .inverseProduct(
// (m_matrix.lowerWithUnitDiag() * (m_matrix.diagonal().asDiagonal())).lower().inverseProduct(vecB));
}

7
Eigen/src/Core/Functors.h
View File

@ -158,7 +158,8 @@ struct ei_functor_traits<ei_scalar_opposite_op<Scalar> >
* \sa class CwiseUnaryOp, MatrixBase::cwiseAbs
*/
template<typename Scalar> struct ei_scalar_abs_op EIGEN_EMPTY_STRUCT {
const Scalar operator() (const Scalar& a) const { return ei_abs(a); }
typedef typename NumTraits<Scalar>::Real result_type;
const result_type operator() (const Scalar& a) const { return ei_abs(a); }
};
template<typename Scalar>
struct ei_functor_traits<ei_scalar_abs_op<Scalar> >
@ -170,8 +171,8 @@ struct ei_functor_traits<ei_scalar_abs_op<Scalar> >
* \sa class CwiseUnaryOp, MatrixBase::cwiseAbs2
*/
template<typename Scalar> struct ei_scalar_abs2_op EIGEN_EMPTY_STRUCT {
const Scalar operator() (const Scalar& a) const { return ei_abs2(a); }
enum { Cost = NumTraits<Scalar>::MulCost };
typedef typename NumTraits<Scalar>::Real result_type;
const result_type operator() (const Scalar& a) const { return ei_abs2(a); }
};
template<typename Scalar>
struct ei_functor_traits<ei_scalar_abs2_op<Scalar> >

37
Eigen/src/Core/MathFunctions.h
View File

@ -181,6 +181,43 @@ inline std::complex<double> ei_exp(std::complex<double> x) { return std::exp(x)
inline std::complex<double> ei_sin(std::complex<double> x) { return std::sin(x); }
inline std::complex<double> ei_cos(std::complex<double> x) { return std::cos(x); }
template<typename T>
inline std::complex<T> ei_sqrt(const std::complex<T>& x)
{
if (std::real(x) == 0.0 && std::imag(x) == 0.0)
return std::complex<T>(0);
else
{
T a = ei_abs(std::real(x));
T b = ei_abs(std::imag(x));
T c;
if (a >= b)
{
T t = b / a;
c = ei_sqrt(a) * ei_sqrt(0.5 * (1.0 + ei_sqrt(1.0 + t * t)));
}
else
{
T t = a / b;
c = ei_sqrt(b) * ei_sqrt(0.5 * (t + ei_sqrt (1.0 + t * t)));
}
T d = std::imag(x) / (2.0 * c);
if (std::real(x) >= 0.0)
{
return std::complex<T>(c, d);
}
else
{
std::complex<T> res(d, c);
if (std::imag(x)<0.0)
res = -res;
return res;
}
}
}
template<> inline std::complex<double> ei_random()
{
return std::complex<double>(ei_random<double>(), ei_random<double>());

22
test/CMakeLists.txt
View File

@ -9,17 +9,17 @@ INCLUDE_DIRECTORIES( ${QT_INCLUDE_DIR} )
SET(test_SRCS
cholesky.cpp
main.cpp
# basicstuff.cpp
# linearstructure.cpp
# product.cpp
# adjoint.cpp
# submatrices.cpp
# miscmatrices.cpp
# smallvectors.cpp
# map.cpp
# cwiseop.cpp
# determinant.cpp
# triangular.cpp
basicstuff.cpp
linearstructure.cpp
product.cpp
adjoint.cpp
submatrices.cpp
miscmatrices.cpp
smallvectors.cpp
map.cpp
cwiseop.cpp
determinant.cpp
triangular.cpp
)
QT4_AUTOMOC(${test_SRCS})

9
test/cholesky.cpp
View File

@ -42,14 +42,15 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
MatrixType a = MatrixType::random(rows,cols).transpose();
VectorType b = VectorType::random(cols);
SquareMatrixType covMat = a.transpose() * a;
SquareMatrixType covMat = a.adjoint() * a;
CholeskyWithoutSquareRoot<SquareMatrixType> cholnosqrt(covMat);
VERIFY_IS_APPROX(covMat, cholnosqrt.matrixU().transpose() * cholnosqrt.vectorD().asDiagonal() * cholnosqrt.matrixU());
VERIFY_IS_APPROX(covMat, cholnosqrt.matrixL() * cholnosqrt.vectorD().asDiagonal() * cholnosqrt.matrixL().adjoint());
VERIFY_IS_APPROX(covMat * cholnosqrt.solve(b), b);
Cholesky<SquareMatrixType> chol(covMat);
VERIFY_IS_APPROX(covMat, chol.matrixU().transpose() * chol.matrixU());
VERIFY_IS_APPROX(covMat, chol.matrixL() * chol.matrixL().adjoint());
VERIFY_IS_APPROX(covMat * chol.solve(b), b);
}
@ -58,7 +59,7 @@ void EigenTest::testCholesky()
for(int i = 0; i < 1; i++) {
cholesky(Matrix3f());
cholesky(Matrix4d());
cholesky(MatrixXd(17,17));
cholesky(MatrixXcd(7,7));
}
}

22
test/main.h
View File

@ -204,17 +204,17 @@ class EigenTest : public QObject
EigenTest(int repeat) : m_repeat(repeat) {}
private slots:
// void testBasicStuff();
// void testLinearStructure();
// void testProduct();
// void testAdjoint();
// void testSubmatrices();
// void testMiscMatrices();
// void testSmallVectors();
// void testMap();
// void testCwiseops();
// void testDeterminant();
// void testTriangular();
void testBasicStuff();
void testLinearStructure();
void testProduct();
void testAdjoint();
void testSubmatrices();
void testMiscMatrices();
void testSmallVectors();
void testMap();
void testCwiseops();
void testDeterminant();
void testTriangular();
void testCholesky();
protected:
int m_repeat;