devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add QR and Cholesky module instantiations in the lib.

Browse Source 
To try it with the unit tests set the cmake variable TEST_LIB to ON.
This commit is contained in:
Gael Guennebaud
2008-06-14 13:02:41 +00:00
parent 53289a8b64
commit f07f907810
15 changed files with 205 additions and 71 deletions
Download Patch File Download Diff File Expand all files Collapse all files

35
Eigen/src/Cholesky/CholeskyInstantiations.cpp Normal file
View File

@@ -0,0 +1,35 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra. Eigen itself is part of the KDE project.
//
// Copyright (C) 2008 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
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifndef EIGEN_EXTERN_INSTANTIATIONS
#define EIGEN_EXTERN_INSTANTIATIONS
#endif
#include "../../Core"
#undef EIGEN_EXTERN_INSTANTIATIONS
#include "../../Cholesky"
namespace Eigen {
EIGEN_CHOLESKY_MODULE_INSTANTIATE();
}

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

@@ -72,7 +72,8 @@ template<typename MatrixType> class CholeskyWithoutSquareRoot
/** \returns whether the matrix is positive definite */
bool isPositiveDefinite(void) const
{
return m_matrix.diagonal().minCoeff() > Scalar(0);
// FIXME is it really correct ?
return m_matrix.diagonal().real().minCoeff() > RealScalar(0);
}
template<typename Derived>

16
Eigen/src/Core/CoreInstanciations.cpp → Eigen/src/Core/CoreInstantiations.cpp
View File

@@ -22,8 +22,8 @@
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifdef EIGEN_EXTERN_INSTANCIATIONS
#undef EIGEN_EXTERN_INSTANCIATIONS
#ifdef EIGEN_EXTERN_INSTANTIATIONS
#undef EIGEN_EXTERN_INSTANTIATIONS
#endif
#include "../../Core"
@@ -31,17 +31,17 @@
namespace Eigen
{
#define EIGEN_INSTANCIATE_PRODUCT(TYPE) \
#define EIGEN_INSTANTIATE_PRODUCT(TYPE) \
template static void ei_cache_friendly_product<TYPE>( \
int _rows, int _cols, int depth, \
bool _lhsRowMajor, const TYPE* _lhs, int _lhsStride, \
bool _rhsRowMajor, const TYPE* _rhs, int _rhsStride, \
bool resRowMajor, TYPE* res, int resStride)
EIGEN_INSTANCIATE_PRODUCT(float);
EIGEN_INSTANCIATE_PRODUCT(double);
EIGEN_INSTANCIATE_PRODUCT(int);
EIGEN_INSTANCIATE_PRODUCT(std::complex<float>);
EIGEN_INSTANCIATE_PRODUCT(std::complex<double>);
EIGEN_INSTANTIATE_PRODUCT(float);
EIGEN_INSTANTIATE_PRODUCT(double);
EIGEN_INSTANTIATE_PRODUCT(int);
EIGEN_INSTANTIATE_PRODUCT(std::complex<float>);
EIGEN_INSTANTIATE_PRODUCT(std::complex<double>);
}

2
Eigen/src/Core/util/ForwardDeclarations.h
View File

@@ -87,7 +87,7 @@ template<typename Scalar> struct ei_scalar_max_op;
template<typename Scalar> struct ei_scalar_random_op;
template<typename Scalar>
static void ei_cache_friendly_product(
void ei_cache_friendly_product(
int _rows, int _cols, int depth,
bool _lhsRowMajor, const Scalar* _lhs, int _lhsStride,
bool _rhsRowMajor, const Scalar* _rhs, int _rhsStride,

15
Eigen/src/QR/HessenbergDecomposition.h
View File

@@ -60,11 +60,11 @@ template<typename _MatrixType> class HessenbergDecomposition
NestByValue<Block<
MatrixType,SizeMinusOne,SizeMinusOne> > > >::RealReturnType SubDiagonalReturnType;
HessenbergDecomposition()
{}
HessenbergDecomposition(int rows, int cols)
: m_matrix(rows,cols), m_hCoeffs(rows-1)
/** 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)
{}
HessenbergDecomposition(const MatrixType& matrix)
@@ -121,6 +121,7 @@ template<typename _MatrixType> class HessenbergDecomposition
CoeffVectorType m_hCoeffs;
};
#ifndef EIGEN_HIDE_HEAVY_CODE
/** \internal
* Performs a tridiagonal decomposition of \a matA in place.
@@ -223,6 +224,8 @@ HessenbergDecomposition<MatrixType>::matrixQ(void) const
return matQ;
}
#endif // EIGEN_HIDE_HEAVY_CODE
/** constructs and returns the matrix H.
* Note that the matrix H is equivalent to the upper part of the packed matrix
* (including the lower sub-diagonal). Therefore, it might be often sufficient
@@ -233,7 +236,7 @@ typename HessenbergDecomposition<MatrixType>::MatrixType
HessenbergDecomposition<MatrixType>::matrixH(void) const
{
// FIXME should this function (and other similar) rather take a matrix as argument
// and fill it (avoids temporaries)
// and fill it (to avoid temporaries)
int n = m_matrix.rows();
MatrixType matH = m_matrix;
matH.corner(BottomLeft,n-2, n-2).template part<Lower>().setZero();

4
Eigen/src/QR/QR.h
View File

@@ -75,6 +75,8 @@ template<typename MatrixType> class QR
VectorType m_hCoeffs;
};
#ifndef EIGEN_HIDE_HEAVY_CODE
template<typename MatrixType>
void QR<MatrixType>::_compute(const MatrixType& matrix)
{
@@ -157,6 +159,8 @@ MatrixType QR<MatrixType>::matrixQ(void) const
return res;
}
#endif // EIGEN_HIDE_HEAVY_CODE
/** \return the QR decomposition of \c *this.
*
* \sa class QR

8
Eigen/src/QR/QrInstanciations.cpp → Eigen/src/QR/QrInstantiations.cpp
View File

@@ -22,9 +22,11 @@
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
#ifdef EIGEN_EXTERN_INSTANCIATIONS
#undef EIGEN_EXTERN_INSTANCIATIONS
#ifndef EIGEN_EXTERN_INSTANTIATIONS
#define EIGEN_EXTERN_INSTANTIATIONS
#endif
#include "../../Core"
#undef EIGEN_EXTERN_INSTANTIATIONS
#include "../../QR"
@@ -36,4 +38,6 @@ template static void ei_tridiagonal_qr_step(double* , double* , int, int, double
template static void ei_tridiagonal_qr_step(float* , float* , int, int, std::complex<float>* , int);
template static void ei_tridiagonal_qr_step(double* , double* , int, int, std::complex<double>* , int);
EIGEN_QR_MODULE_INSTANTIATE();
}

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

@@ -223,7 +223,7 @@ MatrixBase<Derived>::matrixNorm() const
::matrixNorm(derived());
}
#ifndef EIGEN_EXTERN_INSTANCIATIONS
#ifndef EIGEN_EXTERN_INSTANTIATIONS
template<typename RealScalar, typename Scalar>
static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int start, int end, Scalar* matrixQ, int n)
{

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

@@ -60,11 +60,11 @@ template<typename _MatrixType> class Tridiagonalization
NestByValue<Block<
MatrixType,SizeMinusOne,SizeMinusOne> > > >::RealReturnType SubDiagonalReturnType;
Tridiagonalization()
{}
Tridiagonalization(int rows, int cols)
: m_matrix(rows,cols), m_hCoeffs(rows-1)
/** This constructor initializes a Tridiagonalization object for
* further use with Tridiagonalization::compute()
*/
Tridiagonalization(int size = Size==Dynamic ? 2 : Size)
: m_matrix(size,size), m_hCoeffs(size-1)
{}
Tridiagonalization(const MatrixType& matrix)
@@ -90,7 +90,7 @@ template<typename _MatrixType> class Tridiagonalization
*
* \sa packedMatrix()
*/
CoeffVectorType householderCoefficients(void) const { return m_hCoeffs; }
inline CoeffVectorType householderCoefficients(void) const { return m_hCoeffs; }
/** \returns the internal result of the decomposition.
*
@@ -108,7 +108,7 @@ template<typename _MatrixType> class Tridiagonalization
*
* See LAPACK for further details on this packed storage.
*/
const MatrixType& packedMatrix(void) const { return m_matrix; }
inline const MatrixType& packedMatrix(void) const { return m_matrix; }
MatrixType matrixQ(void) const;
MatrixType matrixT(void) const;
@@ -128,6 +128,44 @@ template<typename _MatrixType> class Tridiagonalization
CoeffVectorType m_hCoeffs;
};
/** \returns an expression of the diagonal vector */
template<typename MatrixType>
const typename Tridiagonalization<MatrixType>::DiagonalReturnType
Tridiagonalization<MatrixType>::diagonal(void) const
{
return m_matrix.diagonal().nestByValue().real();
}
/** \returns an expression of the sub-diagonal vector */
template<typename MatrixType>
const typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
Tridiagonalization<MatrixType>::subDiagonal(void) const
{
int n = m_matrix.rows();
return Block<MatrixType,SizeMinusOne,SizeMinusOne>(m_matrix, 1, 0, n-1,n-1)
.nestByValue().diagonal().nestByValue().real();
}
/** constructs and returns the tridiagonal matrix T.
* Note that the matrix T is equivalent to the diagonal and sub-diagonal of the packed matrix.
* Therefore, it might be often sufficient to directly use the packed matrix, or the vector
* expressions returned by diagonal() and subDiagonal() instead of creating a new matrix.
*/
template<typename MatrixType>
typename Tridiagonalization<MatrixType>::MatrixType
Tridiagonalization<MatrixType>::matrixT(void) const
{
// FIXME should this function (and other similar) rather take a matrix as argument
// and fill it (avoids temporaries)
int n = m_matrix.rows();
MatrixType matT = m_matrix;
matT.corner(TopRight,n-1, n-1).diagonal() = subDiagonal().conjugate();
matT.corner(TopRight,n-2, n-2).template part<Upper>().setZero();
matT.corner(BottomLeft,n-2, n-2).template part<Lower>().setZero();
return matT;
}
#ifndef EIGEN_HIDE_HEAVY_CODE
/** \internal
* Performs a tridiagonal decomposition of \a matA in place.
@@ -235,43 +273,6 @@ Tridiagonalization<MatrixType>::matrixQ(void) const
return matQ;
}
/** \returns an expression of the diagonal vector */
template<typename MatrixType>
const typename Tridiagonalization<MatrixType>::DiagonalReturnType
Tridiagonalization<MatrixType>::diagonal(void) const
{
return m_matrix.diagonal().nestByValue().real();
}
/** \returns an expression of the sub-diagonal vector */
template<typename MatrixType>
const typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
Tridiagonalization<MatrixType>::subDiagonal(void) const
{
int n = m_matrix.rows();
return Block<MatrixType,SizeMinusOne,SizeMinusOne>(m_matrix, 1, 0, n-1,n-1)
.nestByValue().diagonal().nestByValue().real();
}
/** constructs and returns the tridiagonal matrix T.
* Note that the matrix T is equivalent to the diagonal and sub-diagonal of the packed matrix.
* Therefore, it might be often sufficient to directly use the packed matrix, or the vector
* expressions returned by diagonal() and subDiagonal() instead of creating a new matrix.
*/
template<typename MatrixType>
typename Tridiagonalization<MatrixType>::MatrixType
Tridiagonalization<MatrixType>::matrixT(void) const
{
// FIXME should this function (and other similar) rather take a matrix as argument
// and fill it (avoids temporaries)
int n = m_matrix.rows();
MatrixType matT = m_matrix;
matT.corner(TopRight,n-1, n-1).diagonal() = subDiagonal().conjugate();
matT.corner(TopRight,n-2, n-2).template part<Upper>().setZero();
matT.corner(BottomLeft,n-2, n-2).template part<Lower>().setZero();
return matT;
}
/** Performs a full decomposition in place */
template<typename MatrixType>
void Tridiagonalization<MatrixType>::decomposeInPlace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
@@ -337,4 +338,6 @@ void Tridiagonalization<MatrixType>::_decomposeInPlace3x3(MatrixType& mat, Diago
}
}
#endif // EIGEN_HIDE_HEAVY_CODE
#endif // EIGEN_TRIDIAGONALIZATION_H