GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-09-11 08:53:14 +08:00
Code Issues Packages Projects Releases Wiki Activity

compilation fixes

Browse Source
This commit is contained in:
Gael Guennebaud 2009-07-27 13:50:23 +02:00
parent 0590c18555
commit 94cc30180e
8 changed files with 35 additions and 43 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
Eigen/src/Core/products/SelfadjointProduct.h
View File

@ -76,9 +76,6 @@ struct ei_selfadjoint_product<Scalar,MatStorageOrder, ColMajor, AAT, UpLo>
Scalar* blockA = ei_aligned_stack_new(Scalar, kc*mc); Scalar* blockA = ei_aligned_stack_new(Scalar, kc*mc);
Scalar* blockB = ei_aligned_stack_new(Scalar, kc*size*Blocking::PacketSize); Scalar* blockB = ei_aligned_stack_new(Scalar, kc*size*Blocking::PacketSize);
// number of columns which can be processed by packet of nr columns
int packet_cols = (size/Blocking::nr)*Blocking::nr;
// note that the actual rhs is the transpose/adjoint of mat // note that the actual rhs is the transpose/adjoint of mat
typedef ei_conj_helper<NumTraits<Scalar>::IsComplex && !AAT, NumTraits<Scalar>::IsComplex && AAT> Conj; typedef ei_conj_helper<NumTraits<Scalar>::IsComplex && !AAT, NumTraits<Scalar>::IsComplex && AAT> Conj;

18
Eigen/src/LU/LU.h
View File

@ -92,7 +92,7 @@ template<typename MatrixType> class LU
MatrixType::MaxColsAtCompileTime // so it has the same number of rows and at most as many columns. MatrixType::MaxColsAtCompileTime // so it has the same number of rows and at most as many columns.
> ImageResultType; > ImageResultType;
/** /**
* \brief Default Constructor. * \brief Default Constructor.
* *
* The default constructor is useful in cases in which the user intends to * The default constructor is useful in cases in which the user intends to
@ -227,7 +227,7 @@ template<typename MatrixType> class LU
* Example: \include LU_solve.cpp * Example: \include LU_solve.cpp
* Output: \verbinclude LU_solve.out * Output: \verbinclude LU_solve.out
* *
* \sa MatrixBase::solveTriangular(), kernel(), computeKernel(), inverse(), computeInverse() * \sa TriangularView::solve(), kernel(), computeKernel(), inverse(), computeInverse()
*/ */
template<typename OtherDerived, typename ResultType> template<typename OtherDerived, typename ResultType>
bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const; bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
@ -380,7 +380,7 @@ void LU<MatrixType>::compute(const MatrixType& matrix)
const int size = matrix.diagonalSize(); const int size = matrix.diagonalSize();
const int rows = matrix.rows(); const int rows = matrix.rows();
const int cols = matrix.cols(); const int cols = matrix.cols();
// this formula comes from experimenting (see "LU precision tuning" thread on the list) // this formula comes from experimenting (see "LU precision tuning" thread on the list)
// and turns out to be identical to Higham's formula used already in LDLt. // and turns out to be identical to Higham's formula used already in LDLt.
m_precision = machine_epsilon<Scalar>() * size; m_precision = machine_epsilon<Scalar>() * size;
@ -486,8 +486,7 @@ void LU<MatrixType>::computeKernel(KernelMatrixType *result) const
y(-m_lu.corner(TopRight, m_rank, dimker)); y(-m_lu.corner(TopRight, m_rank, dimker));
m_lu.corner(TopLeft, m_rank, m_rank) m_lu.corner(TopLeft, m_rank, m_rank)
.template marked<UpperTriangular>() .template triangularView<UpperTriangular>().solveInPlace(y);
.solveTriangularInPlace(y);
for(int i = 0; i < m_rank; ++i) result->row(m_q.coeff(i)) = y.row(i); for(int i = 0; i < m_rank; ++i) result->row(m_q.coeff(i)) = y.row(i);
for(int i = m_rank; i < cols; ++i) result->row(m_q.coeff(i)).setZero(); for(int i = m_rank; i < cols; ++i) result->row(m_q.coeff(i)).setZero();
@ -552,9 +551,8 @@ bool LU<MatrixType>::solve(
for(int i = 0; i < rows; ++i) c.row(m_p.coeff(i)) = b.row(i); for(int i = 0; i < rows; ++i) c.row(m_p.coeff(i)) = b.row(i);
// Step 2 // Step 2
m_lu.corner(Eigen::TopLeft,smalldim,smalldim).template marked<UnitLowerTriangular>() m_lu.corner(Eigen::TopLeft,smalldim,smalldim).template triangularView<UnitLowerTriangular>()
.solveTriangularInPlace( .solveInPlace(c.corner(Eigen::TopLeft, smalldim, c.cols()));
c.corner(Eigen::TopLeft, smalldim, c.cols()));
if(rows>cols) if(rows>cols)
{ {
c.corner(Eigen::BottomLeft, rows-cols, c.cols()) c.corner(Eigen::BottomLeft, rows-cols, c.cols())
@ -572,8 +570,8 @@ bool LU<MatrixType>::solve(
return false; return false;
} }
m_lu.corner(TopLeft, m_rank, m_rank) m_lu.corner(TopLeft, m_rank, m_rank)
.template marked<UpperTriangular>() .template triangularView<UpperTriangular>()
.solveTriangularInPlace(c.corner(TopLeft, m_rank, c.cols())); .solveInPlace(c.corner(TopLeft, m_rank, c.cols()));
// Step 4 // Step 4
result->resize(m_lu.cols(), b.cols()); result->resize(m_lu.cols(), b.cols());

14
Eigen/src/LU/PartialLU.h
View File

@ -70,7 +70,7 @@ template<typename MatrixType> class PartialLU
MatrixType::MaxRowsAtCompileTime) MatrixType::MaxRowsAtCompileTime)
}; };
/** /**
* \brief Default Constructor. * \brief Default Constructor.
* *
* The default constructor is useful in cases in which the user intends to * The default constructor is useful in cases in which the user intends to
@ -99,7 +99,7 @@ template<typename MatrixType> class PartialLU
{ {
ei_assert(m_isInitialized && "PartialLU is not initialized."); ei_assert(m_isInitialized && "PartialLU is not initialized.");
return m_lu; return m_lu;
} }
/** \returns a vector of integers, whose size is the number of rows of the matrix being decomposed, /** \returns a vector of integers, whose size is the number of rows of the matrix being decomposed,
* representing the P permutation i.e. the permutation of the rows. For its precise meaning, * representing the P permutation i.e. the permutation of the rows. For its precise meaning,
@ -127,7 +127,7 @@ template<typename MatrixType> class PartialLU
* Example: \include PartialLU_solve.cpp * Example: \include PartialLU_solve.cpp
* Output: \verbinclude PartialLU_solve.out * Output: \verbinclude PartialLU_solve.out
* *
* \sa MatrixBase::solveTriangular(), inverse(), computeInverse() * \sa TriangularView::solve(), inverse(), computeInverse()
*/ */
template<typename OtherDerived, typename ResultType> template<typename OtherDerived, typename ResultType>
void solve(const MatrixBase<OtherDerived>& b, ResultType *result) const; void solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
@ -235,7 +235,7 @@ void PartialLU<MatrixType>::compute(const MatrixType& matrix)
* one path * one path
*/ */
for(int col = k + 1; col < size; ++col) for(int col = k + 1; col < size; ++col)
m_lu.col(col).end(size-k-1) -= m_lu.col(k).end(size-k-1) * m_lu.coeff(k,col); m_lu.col(col).end(size-k-1) -= m_lu.col(k).end(size-k-1) * m_lu.coeff(k,col);
} }
} }
@ -280,12 +280,10 @@ void PartialLU<MatrixType>::solve(
for(int i = 0; i < size; ++i) result->row(m_p.coeff(i)) = b.row(i); for(int i = 0; i < size; ++i) result->row(m_p.coeff(i)) = b.row(i);
// Step 2 // Step 2
m_lu.template marked<UnitLowerTriangular>() m_lu.template triangularView<UnitLowerTriangular>().solveInPlace(*result);
.solveTriangularInPlace(*result);
// Step 3 // Step 3
m_lu.template marked<UpperTriangular>() m_lu.template triangularView<UpperTriangular>().solveInPlace(*result);
.solveTriangularInPlace(*result);
} }
/** \lu_module /** \lu_module

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

@ -243,7 +243,7 @@ HessenbergDecomposition<MatrixType>::matrixH(void) const
int n = m_matrix.rows(); int n = m_matrix.rows();
MatrixType matH = m_matrix; MatrixType matH = m_matrix;
if (n>2) if (n>2)
matH.corner(BottomLeft,n-2, n-2).template part<LowerTriangular>().setZero(); matH.corner(BottomLeft,n-2, n-2).template triangularView<LowerTriangular>().setZero();
return matH; return matH;
} }

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

@ -51,7 +51,7 @@ template<typename MatrixType> class HouseholderQR
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> MatrixTypeR; typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> MatrixTypeR;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType; typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
/** /**
* \brief Default Constructor. * \brief Default Constructor.
* *
* The default constructor is useful in cases in which the user intends to * The default constructor is useful in cases in which the user intends to
@ -66,14 +66,14 @@ template<typename MatrixType> class HouseholderQR
{ {
compute(matrix); compute(matrix);
} }
/** \returns a read-only expression of the matrix R of the actual the QR decomposition */ /** \returns a read-only expression of the matrix R of the actual the QR decomposition */
const TriangularView<NestByValue<MatrixRBlockType>, UpperTriangular> const TriangularView<NestByValue<MatrixRBlockType>, UpperTriangular>
matrixR(void) const matrixR(void) const
{ {
ei_assert(m_isInitialized && "HouseholderQR is not initialized."); ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
int cols = m_qr.cols(); int cols = m_qr.cols();
return MatrixRBlockType(m_qr, 0, 0, cols, cols).nestByValue().template part<UpperTriangular>(); return MatrixRBlockType(m_qr, 0, 0, cols, cols).nestByValue().template triangularView<UpperTriangular>();
} }
/** This method finds a solution x to the equation Ax=b, where A is the matrix of which /** This method finds a solution x to the equation Ax=b, where A is the matrix of which
@ -95,7 +95,7 @@ template<typename MatrixType> class HouseholderQR
void solve(const MatrixBase<OtherDerived>& b, ResultType *result) const; void solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
MatrixType matrixQ(void) const; MatrixType matrixQ(void) const;
/** \returns a reference to the matrix where the Householder QR decomposition is stored /** \returns a reference to the matrix where the Householder QR decomposition is stored
* in a LAPACK-compatible way. * in a LAPACK-compatible way.
*/ */
@ -113,7 +113,7 @@ template<typename MatrixType> class HouseholderQR
template<typename MatrixType> template<typename MatrixType>
void HouseholderQR<MatrixType>::compute(const MatrixType& matrix) void HouseholderQR<MatrixType>::compute(const MatrixType& matrix)
{ {
m_qr = matrix; m_qr = matrix;
m_hCoeffs.resize(matrix.cols()); m_hCoeffs.resize(matrix.cols());
@ -185,15 +185,15 @@ void HouseholderQR<MatrixType>::solve(
const int rows = m_qr.rows(); const int rows = m_qr.rows();
ei_assert(b.rows() == rows); ei_assert(b.rows() == rows);
result->resize(rows, b.cols()); result->resize(rows, b.cols());
// TODO(keir): There is almost certainly a faster way to multiply by // TODO(keir): There is almost certainly a faster way to multiply by
// Q^T without explicitly forming matrixQ(). Investigate. // Q^T without explicitly forming matrixQ(). Investigate.
*result = matrixQ().transpose()*b; *result = matrixQ().transpose()*b;
const int rank = std::min(result->rows(), result->cols()); const int rank = std::min(result->rows(), result->cols());
m_qr.corner(TopLeft, rank, rank) m_qr.corner(TopLeft, rank, rank)
.template marked<UpperTriangular>() .template triangularView<UpperTriangular>()
.solveTriangularInPlace(result->corner(TopLeft, rank, result->cols())); .solveInPlace(result->corner(TopLeft, rank, result->cols()));
} }
/** \returns the matrix Q */ /** \returns the matrix Q */

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

@ -173,8 +173,8 @@ Tridiagonalization<MatrixType>::matrixT(void) const
matT.corner(TopRight,n-1, n-1).diagonal() = subDiagonal().template cast<Scalar>().conjugate(); matT.corner(TopRight,n-1, n-1).diagonal() = subDiagonal().template cast<Scalar>().conjugate();
if (n>2) if (n>2)
{ {
matT.corner(TopRight,n-2, n-2).template part<UpperTriangular>().setZero(); matT.corner(TopRight,n-2, n-2).template triangularView<UpperTriangular>().setZero();
matT.corner(BottomLeft,n-2, n-2).template part<LowerTriangular>().setZero(); matT.corner(BottomLeft,n-2, n-2).template triangularView<LowerTriangular>().setZero();
} }
return matT; return matT;
} }

3
test/product_syrk.cpp
View File

@ -42,8 +42,7 @@ template<typename MatrixType> void syrk(const MatrixType& m)
Rhs2 rhs2 = Rhs2::Random(rows, ei_random<int>(1,320)); Rhs2 rhs2 = Rhs2::Random(rows, ei_random<int>(1,320));
Rhs3 rhs3 = Rhs3::Random(ei_random<int>(1,320), rows); Rhs3 rhs3 = Rhs3::Random(ei_random<int>(1,320), rows);
Scalar s1 = ei_random<Scalar>(), Scalar s1 = ei_random<Scalar>();
s2 = ei_random<Scalar>();
m2.setZero(); m2.setZero();
VERIFY_IS_APPROX((m2.template selfadjointView<LowerTriangular>().rankUpdate(rhs2,s1)._expression()), VERIFY_IS_APPROX((m2.template selfadjointView<LowerTriangular>().rankUpdate(rhs2,s1)._expression()),

16
test/product_trmv.cpp
View File

@ -24,7 +24,7 @@
#include "main.h" #include "main.h"
template<typename MatrixType> void product_triangular(const MatrixType& m) template<typename MatrixType> void trmv(const MatrixType& m)
{ {
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar; typedef typename NumTraits<Scalar>::Real RealScalar;
@ -79,14 +79,14 @@ template<typename MatrixType> void product_triangular(const MatrixType& m)
// TODO check with sub-matrices // TODO check with sub-matrices
} }
void test_product_triangular() void test_product_trmv()
{ {
for(int i = 0; i < g_repeat ; i++) { for(int i = 0; i < g_repeat ; i++) {
CALL_SUBTEST( product_triangular(Matrix<float, 1, 1>()) ); CALL_SUBTEST( trmv(Matrix<float, 1, 1>()) );
CALL_SUBTEST( product_triangular(Matrix<float, 2, 2>()) ); CALL_SUBTEST( trmv(Matrix<float, 2, 2>()) );
CALL_SUBTEST( product_triangular(Matrix3d()) ); CALL_SUBTEST( trmv(Matrix3d()) );
CALL_SUBTEST( product_triangular(Matrix<std::complex<float>,23, 23>()) ); CALL_SUBTEST( trmv(Matrix<std::complex<float>,23, 23>()) );
CALL_SUBTEST( product_triangular(MatrixXcd(17,17)) ); CALL_SUBTEST( trmv(MatrixXcd(17,17)) );
CALL_SUBTEST( product_triangular(Matrix<float,Dynamic,Dynamic,RowMajor>(19, 19)) ); CALL_SUBTEST( trmv(Matrix<float,Dynamic,Dynamic,RowMajor>(19, 19)) );
} }
} }