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 template to specify QR permutation index type, Fix ColPivHouseholderQR Lapacke bindings

Browse Source
This commit is contained in:
Charles Schlosser
2023-01-11 15:57:28 +00:00
parent 9463fc95f4
commit be7791e097
10 changed files with 309 additions and 252 deletions
Download Patch File Download Diff File Expand all files Collapse all files

105
test/qr_colpivoting.cpp
View File

@@ -13,9 +13,8 @@
#include <Eigen/SVD>
#include "solverbase.h"
template <typename MatrixType>
template <typename MatrixType, typename StorageIndex>
void cod() {
STATIC_CHECK(( internal::is_same<typename CompleteOrthogonalDecomposition<MatrixType>::StorageIndex,int>::value ));
Index rows = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE);
Index cols = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE);
@@ -28,7 +27,7 @@ void cod() {
MatrixQType;
MatrixType matrix;
createRandomPIMatrixOfRank(rank, rows, cols, matrix);
CompleteOrthogonalDecomposition<MatrixType> cod(matrix);
CompleteOrthogonalDecomposition<MatrixType, StorageIndex> cod(matrix);
VERIFY(rank == cod.rank());
VERIFY(cols - cod.rank() == cod.dimensionOfKernel());
VERIFY(!cod.isInjective());
@@ -63,14 +62,14 @@ void cod() {
VERIFY_IS_APPROX(cod_solution, pinv * rhs);
}
template <typename MatrixType, int Cols2>
template <typename MatrixType, typename StorageIndex, int Cols2>
void cod_fixedsize() {
enum {
Rows = MatrixType::RowsAtCompileTime,
Cols = MatrixType::ColsAtCompileTime
};
typedef typename MatrixType::Scalar Scalar;
typedef CompleteOrthogonalDecomposition<Matrix<Scalar, Rows, Cols> > COD;
typedef CompleteOrthogonalDecomposition<Matrix<Scalar, Rows, Cols>, StorageIndex> COD;
int rank = internal::random<int>(1, (std::min)(int(Rows), int(Cols)) - 1);
Matrix<Scalar, Rows, Cols> matrix;
createRandomPIMatrixOfRank(rank, Rows, Cols, matrix);
@@ -96,12 +95,10 @@ void cod_fixedsize() {
VERIFY_IS_APPROX(cod_solution, pinv * rhs);
}
template<typename MatrixType> void qr()
template<typename MatrixType, typename StorageIndex> void qr()
{
using std::sqrt;
STATIC_CHECK(( internal::is_same<typename ColPivHouseholderQR<MatrixType>::StorageIndex,int>::value ));
Index rows = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE), cols = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE), cols2 = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE);
Index rank = internal::random<Index>(1, (std::min)(rows, cols)-1);
@@ -110,7 +107,7 @@ template<typename MatrixType> void qr()
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType;
MatrixType m1;
createRandomPIMatrixOfRank(rank,rows,cols,m1);
ColPivHouseholderQR<MatrixType> qr(m1);
ColPivHouseholderQR<MatrixType, StorageIndex> qr(m1);
VERIFY_IS_EQUAL(rank, qr.rank());
VERIFY_IS_EQUAL(cols - qr.rank(), qr.dimensionOfKernel());
VERIFY(!qr.isInjective());
@@ -158,7 +155,7 @@ template<typename MatrixType> void qr()
}
}
template<typename MatrixType, int Cols2> void qr_fixedsize()
template<typename MatrixType, typename StorageIndex, int Cols2> void qr_fixedsize()
{
using std::sqrt;
using std::abs;
@@ -168,7 +165,7 @@ template<typename MatrixType, int Cols2> void qr_fixedsize()
int rank = internal::random<int>(1, (std::min)(int(Rows), int(Cols))-1);
Matrix<Scalar,Rows,Cols> m1;
createRandomPIMatrixOfRank(rank,Rows,Cols,m1);
ColPivHouseholderQR<Matrix<Scalar,Rows,Cols> > qr(m1);
ColPivHouseholderQR<Matrix<Scalar,Rows,Cols>, StorageIndex> qr(m1);
VERIFY_IS_EQUAL(rank, qr.rank());
VERIFY_IS_EQUAL(Cols - qr.rank(), qr.dimensionOfKernel());
VERIFY_IS_EQUAL(qr.isInjective(), (rank == Rows));
@@ -207,7 +204,7 @@ template<typename MatrixType, int Cols2> void qr_fixedsize()
// for rank-revealing QR. See
// http://www.netlib.org/lapack/lawnspdf/lawn176.pdf
// page 3 for more detail.
template<typename MatrixType> void qr_kahan_matrix()
template<typename MatrixType, typename StorageIndex> void qr_kahan_matrix()
{
using std::sqrt;
using std::abs;
@@ -227,7 +224,7 @@ template<typename MatrixType> void qr_kahan_matrix()
pow_s_i *= s;
}
m1 = (m1 + m1.transpose()).eval();
ColPivHouseholderQR<MatrixType> qr(m1);
ColPivHouseholderQR<MatrixType, StorageIndex> qr(m1);
MatrixType r = qr.matrixQR().template triangularView<Upper>();
RealScalar threshold =
@@ -247,7 +244,7 @@ template<typename MatrixType> void qr_kahan_matrix()
}
}
template<typename MatrixType> void qr_invertible()
template<typename MatrixType, typename StorageIndex> void qr_invertible()
{
using std::log;
using std::abs;
@@ -266,7 +263,7 @@ template<typename MatrixType> void qr_invertible()
m1 += a * a.adjoint();
}
ColPivHouseholderQR<MatrixType> qr(m1);
ColPivHouseholderQR<MatrixType, StorageIndex> qr(m1);
check_solverbase<MatrixType, MatrixType>(m1, qr, size, size, size);
@@ -283,11 +280,11 @@ template<typename MatrixType> void qr_invertible()
VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant());
}
template<typename MatrixType> void qr_verify_assert()
template<typename MatrixType, typename StorageIndex> void qr_verify_assert()
{
MatrixType tmp;
ColPivHouseholderQR<MatrixType> qr;
ColPivHouseholderQR<MatrixType, StorageIndex> qr;
VERIFY_RAISES_ASSERT(qr.matrixQR())
VERIFY_RAISES_ASSERT(qr.solve(tmp))
VERIFY_RAISES_ASSERT(qr.transpose().solve(tmp))
@@ -303,11 +300,11 @@ template<typename MatrixType> void qr_verify_assert()
VERIFY_RAISES_ASSERT(qr.logAbsDeterminant())
}
template<typename MatrixType> void cod_verify_assert()
template<typename MatrixType, typename StorageIndex> void cod_verify_assert()
{
MatrixType tmp;
CompleteOrthogonalDecomposition<MatrixType> cod;
CompleteOrthogonalDecomposition<MatrixType, StorageIndex> cod;
VERIFY_RAISES_ASSERT(cod.matrixQTZ())
VERIFY_RAISES_ASSERT(cod.solve(tmp))
VERIFY_RAISES_ASSERT(cod.transpose().solve(tmp))
@@ -323,50 +320,58 @@ template<typename MatrixType> void cod_verify_assert()
VERIFY_RAISES_ASSERT(cod.logAbsDeterminant())
}
EIGEN_DECLARE_TEST(qr_colpivoting)
{
#if defined(EIGEN_USE_LAPACKE)
typedef lapack_int StorageIndex;
#else
typedef int StorageIndex;
#endif
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1( qr<MatrixXf>() );
CALL_SUBTEST_2( qr<MatrixXd>() );
CALL_SUBTEST_3( qr<MatrixXcd>() );
CALL_SUBTEST_4(( qr_fixedsize<Matrix<float,3,5>, 4 >() ));
CALL_SUBTEST_5(( qr_fixedsize<Matrix<double,6,2>, 3 >() ));
CALL_SUBTEST_5(( qr_fixedsize<Matrix<double,1,1>, 1 >() ));
CALL_SUBTEST_1( (qr<MatrixXf, StorageIndex>)() );
CALL_SUBTEST_2( (qr<MatrixXf, StorageIndex>)() );
CALL_SUBTEST_3( (qr<MatrixXf, StorageIndex>)() );
CALL_SUBTEST_4(( (qr_fixedsize<Matrix<float,3,5>, StorageIndex, 4>)() ));
CALL_SUBTEST_5(( (qr_fixedsize<Matrix<double,6,2>, StorageIndex, 3>)() ));
CALL_SUBTEST_5(( (qr_fixedsize<Matrix<double,1,1>, StorageIndex, 1>)() ));
}
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1( cod<MatrixXf>() );
CALL_SUBTEST_2( cod<MatrixXd>() );
CALL_SUBTEST_3( cod<MatrixXcd>() );
CALL_SUBTEST_4(( cod_fixedsize<Matrix<float,3,5>, 4 >() ));
CALL_SUBTEST_5(( cod_fixedsize<Matrix<double,6,2>, 3 >() ));
CALL_SUBTEST_5(( cod_fixedsize<Matrix<double,1,1>, 1 >() ));
CALL_SUBTEST_1( (cod<MatrixXf, StorageIndex>)() );
CALL_SUBTEST_2( (cod<MatrixXf, StorageIndex>)() );
CALL_SUBTEST_3( (cod<MatrixXf, StorageIndex>)() );
CALL_SUBTEST_4(( (cod_fixedsize<Matrix<float,3,5>, StorageIndex, 4>)() ));
CALL_SUBTEST_5(( (cod_fixedsize<Matrix<double,6,2>, StorageIndex, 3>)() ));
CALL_SUBTEST_5(( (cod_fixedsize<Matrix<double,1,1>, StorageIndex, 1>)() ));
}
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1( qr_invertible<MatrixXf>() );
CALL_SUBTEST_2( qr_invertible<MatrixXd>() );
CALL_SUBTEST_6( qr_invertible<MatrixXcf>() );
CALL_SUBTEST_3( qr_invertible<MatrixXcd>() );
CALL_SUBTEST_1( (qr_invertible<MatrixXf, StorageIndex>)() );
CALL_SUBTEST_2( (qr_invertible<MatrixXd, StorageIndex>)() );
CALL_SUBTEST_6( (qr_invertible<MatrixXcf, StorageIndex>)() );
CALL_SUBTEST_3( (qr_invertible<MatrixXcd, StorageIndex>)() );
}
CALL_SUBTEST_7(qr_verify_assert<Matrix3f>());
CALL_SUBTEST_8(qr_verify_assert<Matrix3d>());
CALL_SUBTEST_1(qr_verify_assert<MatrixXf>());
CALL_SUBTEST_2(qr_verify_assert<MatrixXd>());
CALL_SUBTEST_6(qr_verify_assert<MatrixXcf>());
CALL_SUBTEST_3(qr_verify_assert<MatrixXcd>());
CALL_SUBTEST_7( (qr_verify_assert<Matrix3f, StorageIndex>) ());
CALL_SUBTEST_8( (qr_verify_assert<Matrix3d, StorageIndex>)());
CALL_SUBTEST_1( (qr_verify_assert<MatrixXf, StorageIndex>)());
CALL_SUBTEST_2( (qr_verify_assert<MatrixXd, StorageIndex>)());
CALL_SUBTEST_6( (qr_verify_assert<MatrixXcf, StorageIndex>)());
CALL_SUBTEST_3( (qr_verify_assert<MatrixXcd, StorageIndex>)());
CALL_SUBTEST_7(cod_verify_assert<Matrix3f>());
CALL_SUBTEST_8(cod_verify_assert<Matrix3d>());
CALL_SUBTEST_1(cod_verify_assert<MatrixXf>());
CALL_SUBTEST_2(cod_verify_assert<MatrixXd>());
CALL_SUBTEST_6(cod_verify_assert<MatrixXcf>());
CALL_SUBTEST_3(cod_verify_assert<MatrixXcd>());
CALL_SUBTEST_7( (cod_verify_assert<Matrix3f, StorageIndex>)());
CALL_SUBTEST_8( (cod_verify_assert<Matrix3d, StorageIndex>)());
CALL_SUBTEST_1( (cod_verify_assert<MatrixXf, StorageIndex>)());
CALL_SUBTEST_2( (cod_verify_assert<MatrixXd, StorageIndex>)());
CALL_SUBTEST_6( (cod_verify_assert<MatrixXcf, StorageIndex>)());
CALL_SUBTEST_3( (cod_verify_assert<MatrixXcd, StorageIndex>)());
// Test problem size constructors
CALL_SUBTEST_9(ColPivHouseholderQR<MatrixXf>(10, 20));
CALL_SUBTEST_9((ColPivHouseholderQR<MatrixXf, StorageIndex>(10, 20)));
CALL_SUBTEST_1( qr_kahan_matrix<MatrixXf>() );
CALL_SUBTEST_2( qr_kahan_matrix<MatrixXd>() );
CALL_SUBTEST_1( (qr_kahan_matrix<MatrixXf, StorageIndex>)() );
CALL_SUBTEST_2( (qr_kahan_matrix<MatrixXd, StorageIndex>)() );
}

57
test/qr_fullpivoting.cpp
View File

@@ -12,9 +12,8 @@
#include <Eigen/QR>
#include "solverbase.h"
template<typename MatrixType> void qr()
template<typename MatrixType, typename StorageIndex> void qr()
{
STATIC_CHECK(( internal::is_same<typename FullPivHouseholderQR<MatrixType>::StorageIndex,int>::value ));
static const int Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime;
Index max_size = EIGEN_TEST_MAX_SIZE;
@@ -28,7 +27,7 @@ template<typename MatrixType> void qr()
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType;
MatrixType m1;
createRandomPIMatrixOfRank(rank,rows,cols,m1);
FullPivHouseholderQR<MatrixType> qr(m1);
FullPivHouseholderQR<MatrixType, StorageIndex> qr(m1);
VERIFY_IS_EQUAL(rank, qr.rank());
VERIFY_IS_EQUAL(cols - qr.rank(), qr.dimensionOfKernel());
VERIFY(!qr.isInjective());
@@ -67,7 +66,7 @@ template<typename MatrixType> void qr()
}
}
template<typename MatrixType> void qr_invertible()
template<typename MatrixType, typename StorageIndex> void qr_invertible()
{
using std::log;
using std::abs;
@@ -88,7 +87,7 @@ template<typename MatrixType> void qr_invertible()
m1 += a * a.adjoint();
}
FullPivHouseholderQR<MatrixType> qr(m1);
FullPivHouseholderQR<MatrixType, StorageIndex> qr(m1);
VERIFY(qr.isInjective());
VERIFY(qr.isInvertible());
VERIFY(qr.isSurjective());
@@ -108,11 +107,11 @@ template<typename MatrixType> void qr_invertible()
VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant());
}
template<typename MatrixType> void qr_verify_assert()
template<typename MatrixType, typename StorageIndex> void qr_verify_assert()
{
MatrixType tmp;
FullPivHouseholderQR<MatrixType> qr;
FullPivHouseholderQR<MatrixType, StorageIndex> qr;
VERIFY_RAISES_ASSERT(qr.matrixQR())
VERIFY_RAISES_ASSERT(qr.solve(tmp))
VERIFY_RAISES_ASSERT(qr.transpose().solve(tmp))
@@ -130,33 +129,35 @@ template<typename MatrixType> void qr_verify_assert()
EIGEN_DECLARE_TEST(qr_fullpivoting)
{
typedef int StorageIndex;
for(int i = 0; i < 1; i++) {
CALL_SUBTEST_5( qr<Matrix3f>() );
CALL_SUBTEST_6( qr<Matrix3d>() );
CALL_SUBTEST_8( qr<Matrix2f>() );
CALL_SUBTEST_1( qr<MatrixXf>() );
CALL_SUBTEST_2( qr<MatrixXd>() );
CALL_SUBTEST_3( qr<MatrixXcd>() );
CALL_SUBTEST_5( (qr<Matrix3f, StorageIndex>()) );
CALL_SUBTEST_6( (qr<Matrix3d, StorageIndex>()) );
CALL_SUBTEST_8( (qr<Matrix2f, StorageIndex>()) );
CALL_SUBTEST_1( (qr<MatrixXf, StorageIndex>()) );
CALL_SUBTEST_2( (qr<MatrixXd, StorageIndex>()) );
CALL_SUBTEST_3( (qr<MatrixXcd, StorageIndex>()) );
}
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1( qr_invertible<MatrixXf>() );
CALL_SUBTEST_2( qr_invertible<MatrixXd>() );
CALL_SUBTEST_4( qr_invertible<MatrixXcf>() );
CALL_SUBTEST_3( qr_invertible<MatrixXcd>() );
CALL_SUBTEST_1( (qr_invertible<MatrixXf, StorageIndex>()) );
CALL_SUBTEST_2( (qr_invertible<MatrixXd, StorageIndex>()) );
CALL_SUBTEST_4( (qr_invertible<MatrixXcf, StorageIndex>()) );
CALL_SUBTEST_3( (qr_invertible<MatrixXcd, StorageIndex>()) );
}
CALL_SUBTEST_5(qr_verify_assert<Matrix3f>());
CALL_SUBTEST_6(qr_verify_assert<Matrix3d>());
CALL_SUBTEST_1(qr_verify_assert<MatrixXf>());
CALL_SUBTEST_2(qr_verify_assert<MatrixXd>());
CALL_SUBTEST_4(qr_verify_assert<MatrixXcf>());
CALL_SUBTEST_3(qr_verify_assert<MatrixXcd>());
CALL_SUBTEST_5( (qr_verify_assert<Matrix3f, StorageIndex>()) );
CALL_SUBTEST_6( (qr_verify_assert<Matrix3d, StorageIndex>()) );
CALL_SUBTEST_1( (qr_verify_assert<MatrixXf, StorageIndex>()) );
CALL_SUBTEST_2( (qr_verify_assert<MatrixXd, StorageIndex>()) );
CALL_SUBTEST_4( (qr_verify_assert<MatrixXcf, StorageIndex>()) );
CALL_SUBTEST_3( (qr_verify_assert<MatrixXcd, StorageIndex>()) );
// Test problem size constructors
CALL_SUBTEST_7(FullPivHouseholderQR<MatrixXf>(10, 20));
CALL_SUBTEST_7((FullPivHouseholderQR<Matrix<float,10,20> >(10,20)));
CALL_SUBTEST_7((FullPivHouseholderQR<Matrix<float,10,20> >(Matrix<float,10,20>::Random())));
CALL_SUBTEST_7((FullPivHouseholderQR<Matrix<float,20,10> >(20,10)));
CALL_SUBTEST_7((FullPivHouseholderQR<Matrix<float,20,10> >(Matrix<float,20,10>::Random())));
CALL_SUBTEST_7( (FullPivHouseholderQR<MatrixXf,StorageIndex>(10, 20)));
CALL_SUBTEST_7( (FullPivHouseholderQR<Matrix<float, 10, 20>, StorageIndex>(10, 20)));
CALL_SUBTEST_7( (FullPivHouseholderQR<Matrix<float, 10, 20>, StorageIndex>(Matrix<float,10,20>::Random())));
CALL_SUBTEST_7( (FullPivHouseholderQR<Matrix<float, 20, 10>, StorageIndex>(20, 10)));
CALL_SUBTEST_7( (FullPivHouseholderQR<Matrix<float, 20, 10>, StorageIndex>(Matrix<float,20,10>::Random())));
}