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

the Index types change.

Browse Source 
As discussed on the list (too long to explain here).
This commit is contained in:
Benoit Jacob
2010-05-30 16:00:58 -04:00
parent faa3ff3be6
commit aaaade4b3d
158 changed files with 3137 additions and 2878 deletions
Download Patch File Download Diff File Expand all files Collapse all files

82
Eigen/src/LU/FullPivLU.h
View File

@@ -68,8 +68,10 @@ template<typename _MatrixType> class FullPivLU
};
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef typename ei_plain_row_type<MatrixType, int>::type IntRowVectorType;
typedef typename ei_plain_col_type<MatrixType, int>::type IntColVectorType;
typedef typename ei_traits<MatrixType>::StorageKind StorageKind;
typedef typename ei_index<StorageKind>::type Index;
typedef typename ei_plain_row_type<MatrixType, Index>::type IntRowVectorType;
typedef typename ei_plain_col_type<MatrixType, Index>::type IntColVectorType;
typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationQType;
typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationPType;
@@ -87,7 +89,7 @@ template<typename _MatrixType> class FullPivLU
* according to the specified problem \a size.
* \sa FullPivLU()
*/
FullPivLU(int rows, int cols);
FullPivLU(Index rows, Index cols);
/** Constructor.
*
@@ -124,7 +126,7 @@ template<typename _MatrixType> class FullPivLU
*
* \sa rank()
*/
inline int nonzeroPivots() const
inline Index nonzeroPivots() const
{
ei_assert(m_isInitialized && "LU is not initialized.");
return m_nonzero_pivots;
@@ -301,12 +303,12 @@ template<typename _MatrixType> class FullPivLU
* For that, it uses the threshold value that you can control by calling
* setThreshold(const RealScalar&).
*/
inline int rank() const
inline Index rank() const
{
ei_assert(m_isInitialized && "LU is not initialized.");
RealScalar premultiplied_threshold = ei_abs(m_maxpivot) * threshold();
int result = 0;
for(int i = 0; i < m_nonzero_pivots; ++i)
Index result = 0;
for(Index i = 0; i < m_nonzero_pivots; ++i)
result += (ei_abs(m_lu.coeff(i,i)) > premultiplied_threshold);
return result;
}
@@ -317,7 +319,7 @@ template<typename _MatrixType> class FullPivLU
* For that, it uses the threshold value that you can control by calling
* setThreshold(const RealScalar&).
*/
inline int dimensionOfKernel() const
inline Index dimensionOfKernel() const
{
ei_assert(m_isInitialized && "LU is not initialized.");
return cols() - rank();
@@ -378,8 +380,8 @@ template<typename _MatrixType> class FullPivLU
MatrixType reconstructedMatrix() const;
inline int rows() const { return m_lu.rows(); }
inline int cols() const { return m_lu.cols(); }
inline Index rows() const { return m_lu.rows(); }
inline Index cols() const { return m_lu.cols(); }
protected:
MatrixType m_lu;
@@ -387,7 +389,7 @@ template<typename _MatrixType> class FullPivLU
PermutationQType m_q;
IntColVectorType m_rowsTranspositions;
IntRowVectorType m_colsTranspositions;
int m_det_pq, m_nonzero_pivots;
Index m_det_pq, m_nonzero_pivots;
RealScalar m_maxpivot, m_prescribedThreshold;
bool m_isInitialized, m_usePrescribedThreshold;
};
@@ -399,7 +401,7 @@ FullPivLU<MatrixType>::FullPivLU()
}
template<typename MatrixType>
FullPivLU<MatrixType>::FullPivLU(int rows, int cols)
FullPivLU<MatrixType>::FullPivLU(Index rows, Index cols)
: m_lu(rows, cols),
m_p(rows),
m_q(cols),
@@ -429,26 +431,26 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
m_isInitialized = true;
m_lu = matrix;
const int size = matrix.diagonalSize();
const int rows = matrix.rows();
const int cols = matrix.cols();
const Index size = matrix.diagonalSize();
const Index rows = matrix.rows();
const Index cols = matrix.cols();
// will store the transpositions, before we accumulate them at the end.
// can't accumulate on-the-fly because that will be done in reverse order for the rows.
m_rowsTranspositions.resize(matrix.rows());
m_colsTranspositions.resize(matrix.cols());
int number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. m_rowsTranspositions[i]!=i
Index number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. m_rowsTranspositions[i]!=i
m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
m_maxpivot = RealScalar(0);
RealScalar cutoff(0);
for(int k = 0; k < size; ++k)
for(Index k = 0; k < size; ++k)
{
// First, we need to find the pivot.
// biggest coefficient in the remaining bottom-right corner (starting at row k, col k)
int row_of_biggest_in_corner, col_of_biggest_in_corner;
Index row_of_biggest_in_corner, col_of_biggest_in_corner;
RealScalar biggest_in_corner;
biggest_in_corner = m_lu.bottomRightCorner(rows-k, cols-k)
.cwiseAbs()
@@ -468,7 +470,7 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
// before exiting, make sure to initialize the still uninitialized transpositions
// in a sane state without destroying what we already have.
m_nonzero_pivots = k;
for(int i = k; i < size; ++i)
for(Index i = k; i < size; ++i)
{
m_rowsTranspositions.coeffRef(i) = i;
m_colsTranspositions.coeffRef(i) = i;
@@ -505,11 +507,11 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
// permutations P and Q
m_p.setIdentity(rows);
for(int k = size-1; k >= 0; --k)
for(Index k = size-1; k >= 0; --k)
m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));
m_q.setIdentity(cols);
for(int k = 0; k < size; ++k)
for(Index k = 0; k < size; ++k)
m_q.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
m_det_pq = (number_of_transpositions%2) ? -1 : 1;
@@ -531,7 +533,7 @@ template<typename MatrixType>
MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const
{
ei_assert(m_isInitialized && "LU is not initialized.");
const int smalldim = std::min(m_lu.rows(), m_lu.cols());
const Index smalldim = std::min(m_lu.rows(), m_lu.cols());
// LU
MatrixType res(m_lu.rows(),m_lu.cols());
// FIXME the .toDenseMatrix() should not be needed...
@@ -564,7 +566,7 @@ struct ei_kernel_retval<FullPivLU<_MatrixType> >
template<typename Dest> void evalTo(Dest& dst) const
{
const int cols = dec().matrixLU().cols(), dimker = cols - rank();
const Index cols = dec().matrixLU().cols(), dimker = cols - rank();
if(dimker == 0)
{
// The Kernel is just {0}, so it doesn't have a basis properly speaking, but let's
@@ -590,10 +592,10 @@ struct ei_kernel_retval<FullPivLU<_MatrixType> >
* independent vectors in Ker U.
*/
Matrix<int, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
int p = 0;
for(int i = 0; i < dec().nonzeroPivots(); ++i)
Index p = 0;
for(Index i = 0; i < dec().nonzeroPivots(); ++i)
if(ei_abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
pivots.coeffRef(p++) = i;
ei_internal_assert(p == rank());
@@ -605,14 +607,14 @@ struct ei_kernel_retval<FullPivLU<_MatrixType> >
Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, MatrixType::Options,
MaxSmallDimAtCompileTime, MatrixType::MaxColsAtCompileTime>
m(dec().matrixLU().block(0, 0, rank(), cols));
for(int i = 0; i < rank(); ++i)
for(Index i = 0; i < rank(); ++i)
{
if(i) m.row(i).head(i).setZero();
m.row(i).tail(cols-i) = dec().matrixLU().row(pivots.coeff(i)).tail(cols-i);
}
m.block(0, 0, rank(), rank());
m.block(0, 0, rank(), rank()).template triangularView<StrictlyLower>().setZero();
for(int i = 0; i < rank(); ++i)
for(Index i = 0; i < rank(); ++i)
m.col(i).swap(m.col(pivots.coeff(i)));
// ok, we have our trapezoid matrix, we can apply the triangular solver.
@@ -624,13 +626,13 @@ struct ei_kernel_retval<FullPivLU<_MatrixType> >
);
// now we must undo the column permutation that we had applied!
for(int i = rank()-1; i >= 0; --i)
for(Index i = rank()-1; i >= 0; --i)
m.col(i).swap(m.col(pivots.coeff(i)));
// see the negative sign in the next line, that's what we were talking about above.
for(int i = 0; i < rank(); ++i) dst.row(dec().permutationQ().indices().coeff(i)) = -m.row(i).tail(dimker);
for(int i = rank(); i < cols; ++i) dst.row(dec().permutationQ().indices().coeff(i)).setZero();
for(int k = 0; k < dimker; ++k) dst.coeffRef(dec().permutationQ().indices().coeff(rank()+k), k) = Scalar(1);
for(Index i = 0; i < rank(); ++i) dst.row(dec().permutationQ().indices().coeff(i)) = -m.row(i).tail(dimker);
for(Index i = rank(); i < cols; ++i) dst.row(dec().permutationQ().indices().coeff(i)).setZero();
for(Index k = 0; k < dimker; ++k) dst.coeffRef(dec().permutationQ().indices().coeff(rank()+k), k) = Scalar(1);
}
};
@@ -658,15 +660,15 @@ struct ei_image_retval<FullPivLU<_MatrixType> >
return;
}
Matrix<int, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
int p = 0;
for(int i = 0; i < dec().nonzeroPivots(); ++i)
Index p = 0;
for(Index i = 0; i < dec().nonzeroPivots(); ++i)
if(ei_abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
pivots.coeffRef(p++) = i;
ei_internal_assert(p == rank());
for(int i = 0; i < rank(); ++i)
for(Index i = 0; i < rank(); ++i)
dst.col(i) = originalMatrix().col(dec().permutationQ().indices().coeff(pivots.coeff(i)));
}
};
@@ -689,10 +691,10 @@ struct ei_solve_retval<FullPivLU<_MatrixType>, Rhs>
* Step 4: result = Q * c;
*/
const int rows = dec().rows(), cols = dec().cols(),
const Index rows = dec().rows(), cols = dec().cols(),
nonzero_pivots = dec().nonzeroPivots();
ei_assert(rhs().rows() == rows);
const int smalldim = std::min(rows, cols);
const Index smalldim = std::min(rows, cols);
if(nonzero_pivots == 0)
{
@@ -724,9 +726,9 @@ struct ei_solve_retval<FullPivLU<_MatrixType>, Rhs>
.solveInPlace(c.topRows(nonzero_pivots));
// Step 4
for(int i = 0; i < nonzero_pivots; ++i)
for(Index i = 0; i < nonzero_pivots; ++i)
dst.row(dec().permutationQ().indices().coeff(i)) = c.row(i);
for(int i = nonzero_pivots; i < dec().matrixLU().cols(); ++i)
for(Index i = nonzero_pivots; i < dec().matrixLU().cols(); ++i)
dst.row(dec().permutationQ().indices().coeff(i)).setZero();
}
};

7
Eigen/src/LU/Inverse.h
View File

@@ -281,7 +281,8 @@ struct ei_traits<ei_inverse_impl<MatrixType> >
template<typename MatrixType>
struct ei_inverse_impl : public ReturnByValue<ei_inverse_impl<MatrixType> >
{
typedef typename MatrixType::Nested MatrixTypeNested;
typedef typename MatrixType::Index Index;
typedef typename ei_eval<MatrixType>::type MatrixTypeNested;
typedef typename ei_cleantype<MatrixTypeNested>::type MatrixTypeNestedCleaned;
const MatrixTypeNested m_matrix;
@@ -290,8 +291,8 @@ struct ei_inverse_impl : public ReturnByValue<ei_inverse_impl<MatrixType> >
: m_matrix(matrix)
{}
inline int rows() const { return m_matrix.rows(); }
inline int cols() const { return m_matrix.cols(); }
inline Index rows() const { return m_matrix.rows(); }
inline Index cols() const { return m_matrix.cols(); }
template<typename Dest> inline void evalTo(Dest& dst) const
{

65
Eigen/src/LU/PartialPivLU.h
View File

@@ -71,7 +71,9 @@ template<typename _MatrixType> class PartialPivLU
};
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef typename ei_plain_col_type<MatrixType, int>::type PermutationVectorType;
typedef typename ei_traits<MatrixType>::StorageKind StorageKind;
typedef typename ei_index<StorageKind>::type Index;
typedef typename ei_plain_col_type<MatrixType, Index>::type PermutationVectorType;
typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
@@ -89,7 +91,7 @@ template<typename _MatrixType> class PartialPivLU
* according to the specified problem \a size.
* \sa PartialPivLU()
*/
PartialPivLU(int size);
PartialPivLU(Index size);
/** Constructor.
*
@@ -178,14 +180,14 @@ template<typename _MatrixType> class PartialPivLU
MatrixType reconstructedMatrix() const;
inline int rows() const { return m_lu.rows(); }
inline int cols() const { return m_lu.cols(); }
inline Index rows() const { return m_lu.rows(); }
inline Index cols() const { return m_lu.cols(); }
protected:
MatrixType m_lu;
PermutationType m_p;
PermutationVectorType m_rowsTranspositions;
int m_det_p;
Index m_det_p;
bool m_isInitialized;
};
@@ -200,7 +202,7 @@ PartialPivLU<MatrixType>::PartialPivLU()
}
template<typename MatrixType>
PartialPivLU<MatrixType>::PartialPivLU(int size)
PartialPivLU<MatrixType>::PartialPivLU(Index size)
: m_lu(size, size),
m_p(size),
m_rowsTranspositions(size),
@@ -233,6 +235,7 @@ struct ei_partial_lu_impl
typedef Block<MapLU, Dynamic, Dynamic> MatrixType;
typedef Block<MatrixType,Dynamic,Dynamic> BlockType;
typedef typename MatrixType::RealScalar RealScalar;
typedef typename MatrixType::Index Index;
/** \internal performs the LU decomposition in-place of the matrix \a lu
* using an unblocked algorithm.
@@ -246,14 +249,14 @@ struct ei_partial_lu_impl
* undefined coefficients (to avoid generating inf/nan values). Returns true
* otherwise.
*/
static bool unblocked_lu(MatrixType& lu, int* row_transpositions, int& nb_transpositions)
static bool unblocked_lu(MatrixType& lu, Index* row_transpositions, Index& nb_transpositions)
{
const int rows = lu.rows();
const int size = std::min(lu.rows(),lu.cols());
const Index rows = lu.rows();
const Index size = std::min(lu.rows(),lu.cols());
nb_transpositions = 0;
for(int k = 0; k < size; ++k)
for(Index k = 0; k < size; ++k)
{
int row_of_biggest_in_col;
Index row_of_biggest_in_col;
RealScalar biggest_in_corner
= lu.col(k).tail(rows-k).cwiseAbs().maxCoeff(&row_of_biggest_in_col);
row_of_biggest_in_col += k;
@@ -265,7 +268,7 @@ struct ei_partial_lu_impl
// the blocked_lu code can't guarantee the same.
// before exiting, make sure to initialize the still uninitialized row_transpositions
// in a sane state without destroying what we already have.
for(int i = k; i < size; i++)
for(Index i = k; i < size; i++)
row_transpositions[i] = i;
return false;
}
@@ -280,8 +283,8 @@ struct ei_partial_lu_impl
if(k<rows-1)
{
int rrows = rows-k-1;
int rsize = size-k-1;
Index rrows = rows-k-1;
Index rsize = size-k-1;
lu.col(k).tail(rrows) /= lu.coeff(k,k);
lu.bottomRightCorner(rrows,rsize).noalias() -= lu.col(k).tail(rrows) * lu.row(k).tail(rsize);
}
@@ -306,12 +309,12 @@ struct ei_partial_lu_impl
* 1 - reduce the number of instanciations to the strict minimum
* 2 - avoid infinite recursion of the instanciations with Block<Block<Block<...> > >
*/
static bool blocked_lu(int rows, int cols, Scalar* lu_data, int luStride, int* row_transpositions, int& nb_transpositions, int maxBlockSize=256)
static bool blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, Index* row_transpositions, Index& nb_transpositions, Index maxBlockSize=256)
{
MapLU lu1(lu_data,StorageOrder==RowMajor?rows:luStride,StorageOrder==RowMajor?luStride:cols);
MatrixType lu(lu1,0,0,rows,cols);
const int size = std::min(rows,cols);
const Index size = std::min(rows,cols);
// if the matrix is too small, no blocking:
if(size<=16)
@@ -321,19 +324,19 @@ struct ei_partial_lu_impl
// automatically adjust the number of subdivisions to the size
// of the matrix so that there is enough sub blocks:
int blockSize;
Index blockSize;
{
blockSize = size/8;
blockSize = (blockSize/16)*16;
blockSize = std::min(std::max(blockSize,8), maxBlockSize);
blockSize = std::min(std::max(blockSize,Index(8)), maxBlockSize);
}
nb_transpositions = 0;
for(int k = 0; k < size; k+=blockSize)
for(Index k = 0; k < size; k+=blockSize)
{
int bs = std::min(size-k,blockSize); // actual size of the block
int trows = rows - k - bs; // trailing rows
int tsize = size - k - bs; // trailing size
Index bs = std::min(size-k,blockSize); // actual size of the block
Index trows = rows - k - bs; // trailing rows
Index tsize = size - k - bs; // trailing size
// partition the matrix:
// A00 | A01 | A02
@@ -346,7 +349,7 @@ struct ei_partial_lu_impl
BlockType A21(lu,k+bs,k,trows,bs);
BlockType A22(lu,k+bs,k+bs,trows,tsize);
int nb_transpositions_in_panel;
Index nb_transpositions_in_panel;
// recursively calls the blocked LU algorithm with a very small
// blocking size:
if(!blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
@@ -355,23 +358,23 @@ struct ei_partial_lu_impl
// end quickly with undefined coefficients, just avoid generating inf/nan values.
// before exiting, make sure to initialize the still uninitialized row_transpositions
// in a sane state without destroying what we already have.
for(int i=k; i<size; ++i)
for(Index i=k; i<size; ++i)
row_transpositions[i] = i;
return false;
}
nb_transpositions += nb_transpositions_in_panel;
// update permutations and apply them to A10
for(int i=k; i<k+bs; ++i)
for(Index i=k; i<k+bs; ++i)
{
int piv = (row_transpositions[i] += k);
Index piv = (row_transpositions[i] += k);
A_0.row(i).swap(A_0.row(piv));
}
if(trows)
{
// apply permutations to A_2
for(int i=k;i<k+bs; ++i)
for(Index i=k;i<k+bs; ++i)
A_2.row(i).swap(A_2.row(row_transpositions[i]));
// A12 = A11^-1 A12
@@ -387,7 +390,7 @@ struct ei_partial_lu_impl
/** \internal performs the LU decomposition with partial pivoting in-place.
*/
template<typename MatrixType, typename IntVector>
void ei_partial_lu_inplace(MatrixType& lu, IntVector& row_transpositions, int& nb_transpositions)
void ei_partial_lu_inplace(MatrixType& lu, IntVector& row_transpositions, typename MatrixType::Index& nb_transpositions)
{
ei_assert(lu.cols() == row_transpositions.size());
ei_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1);
@@ -403,16 +406,16 @@ PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const MatrixType& ma
m_lu = matrix;
ei_assert(matrix.rows() == matrix.cols() && "PartialPivLU is only for square (and moreover invertible) matrices");
const int size = matrix.rows();
const Index size = matrix.rows();
m_rowsTranspositions.resize(size);
int nb_transpositions;
Index nb_transpositions;
ei_partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions);
m_det_p = (nb_transpositions%2) ? -1 : 1;
m_p.setIdentity(size);
for(int k = size-1; k >= 0; --k)
for(Index k = size-1; k >= 0; --k)
m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));
m_isInitialized = true;