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

Bug #1796: Make matrix squareroot usable for Map and Ref types

Browse Source
This commit is contained in:
Christoph Hertzberg 2019-12-20 18:10:22 +01:00
parent d86544d654
commit 1e9664b147
2 changed files with 48 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h
View File

@ -253,18 +253,19 @@ struct matrix_sqrt_compute
template <typename MatrixType> template <typename MatrixType>
struct matrix_sqrt_compute<MatrixType, 0> struct matrix_sqrt_compute<MatrixType, 0>
{ {
typedef typename MatrixType::PlainObject PlainType;
template <typename ResultType> template <typename ResultType>
static void run(const MatrixType &arg, ResultType &result) static void run(const MatrixType &arg, ResultType &result)
{ {
eigen_assert(arg.rows() == arg.cols()); eigen_assert(arg.rows() == arg.cols());
// Compute Schur decomposition of arg // Compute Schur decomposition of arg
const RealSchur<MatrixType> schurOfA(arg); const RealSchur<PlainType> schurOfA(arg);
const MatrixType& T = schurOfA.matrixT(); const PlainType& T = schurOfA.matrixT();
const MatrixType& U = schurOfA.matrixU(); const PlainType& U = schurOfA.matrixU();
// Compute square root of T // Compute square root of T
MatrixType sqrtT = MatrixType::Zero(arg.rows(), arg.cols()); PlainType sqrtT = PlainType::Zero(arg.rows(), arg.cols());
matrix_sqrt_quasi_triangular(T, sqrtT); matrix_sqrt_quasi_triangular(T, sqrtT);
// Compute square root of arg // Compute square root of arg
@ -278,18 +279,19 @@ struct matrix_sqrt_compute<MatrixType, 0>
template <typename MatrixType> template <typename MatrixType>
struct matrix_sqrt_compute<MatrixType, 1> struct matrix_sqrt_compute<MatrixType, 1>
{ {
typedef typename MatrixType::PlainObject PlainType;
template <typename ResultType> template <typename ResultType>
static void run(const MatrixType &arg, ResultType &result) static void run(const MatrixType &arg, ResultType &result)
{ {
eigen_assert(arg.rows() == arg.cols()); eigen_assert(arg.rows() == arg.cols());
// Compute Schur decomposition of arg // Compute Schur decomposition of arg
const ComplexSchur<MatrixType> schurOfA(arg); const ComplexSchur<PlainType> schurOfA(arg);
const MatrixType& T = schurOfA.matrixT(); const PlainType& T = schurOfA.matrixT();
const MatrixType& U = schurOfA.matrixU(); const PlainType& U = schurOfA.matrixU();
// Compute square root of T // Compute square root of T
MatrixType sqrtT; PlainType sqrtT;
matrix_sqrt_triangular(T, sqrtT); matrix_sqrt_triangular(T, sqrtT);
// Compute square root of arg // Compute square root of arg

38
unsupported/test/matrix_function.cpp
View File

@ -177,6 +177,39 @@ void testMatrixType(const MatrixType& m)
} }
} }
template<typename MatrixType>
void testMapRef(const MatrixType& A)
{
// Test if passing Ref and Map objects is possible
// (Regression test for Bug #1796)
Index size = A.rows();
MatrixType X; X.setRandom(size, size);
MatrixType Y(size,size);
Ref< MatrixType> R(Y);
Ref<const MatrixType> Rc(X);
Map< MatrixType> M(Y.data(), size, size);
Map<const MatrixType> Mc(X.data(), size, size);
X = X*X; // make sure sqrt is possible
Y = X.sqrt();
R = Rc.sqrt();
M = Mc.sqrt();
Y = X.exp();
R = Rc.exp();
M = Mc.exp();
X = Y; // make sure log is possible
Y = X.log();
R = Rc.log();
M = Mc.log();
Y = X.cos() + Rc.cos() + Mc.cos();
Y = X.sin() + Rc.sin() + Mc.sin();
Y = X.cosh() + Rc.cosh() + Mc.cosh();
Y = X.sinh() + Rc.sinh() + Mc.sinh();
}
EIGEN_DECLARE_TEST(matrix_function) EIGEN_DECLARE_TEST(matrix_function)
{ {
CALL_SUBTEST_1(testMatrixType(Matrix<float,1,1>())); CALL_SUBTEST_1(testMatrixType(Matrix<float,1,1>()));
@ -186,4 +219,9 @@ EIGEN_DECLARE_TEST(matrix_function)
CALL_SUBTEST_5(testMatrixType(Matrix<double,5,5,RowMajor>())); CALL_SUBTEST_5(testMatrixType(Matrix<double,5,5,RowMajor>()));
CALL_SUBTEST_6(testMatrixType(Matrix4cd())); CALL_SUBTEST_6(testMatrixType(Matrix4cd()));
CALL_SUBTEST_7(testMatrixType(MatrixXd(13,13))); CALL_SUBTEST_7(testMatrixType(MatrixXd(13,13)));
CALL_SUBTEST_1(testMapRef(Matrix<float,1,1>()));
CALL_SUBTEST_2(testMapRef(Matrix3cf()));
CALL_SUBTEST_3(testMapRef(MatrixXf(8,8)));
CALL_SUBTEST_7(testMapRef(MatrixXd(13,13)));
} }