diff options
Diffstat (limited to 'unsupported/Eigen/src/MatrixFunctions')
6 files changed, 44 insertions, 52 deletions
diff --git a/unsupported/Eigen/src/MatrixFunctions/CMakeLists.txt b/unsupported/Eigen/src/MatrixFunctions/CMakeLists.txt deleted file mode 100644 index cdde64d2c..000000000 --- a/unsupported/Eigen/src/MatrixFunctions/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_MatrixFunctions_SRCS "*.h") - -INSTALL(FILES - ${Eigen_MatrixFunctions_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/MatrixFunctions COMPONENT Devel - ) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h index bbb7e5776..4bb1852b6 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h @@ -65,7 +65,7 @@ template <typename MatrixType> void matrix_exp_pade3(const MatrixType &A, MatrixType &U, MatrixType &V) { typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar; - const RealScalar b[] = {120., 60., 12., 1.}; + const RealScalar b[] = {120.L, 60.L, 12.L, 1.L}; const MatrixType A2 = A * A; const MatrixType tmp = b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols()); U.noalias() = A * tmp; @@ -81,7 +81,7 @@ template <typename MatrixType> void matrix_exp_pade5(const MatrixType &A, MatrixType &U, MatrixType &V) { typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar; - const RealScalar b[] = {30240., 15120., 3360., 420., 30., 1.}; + const RealScalar b[] = {30240.L, 15120.L, 3360.L, 420.L, 30.L, 1.L}; const MatrixType A2 = A * A; const MatrixType A4 = A2 * A2; const MatrixType tmp = b[5] * A4 + b[3] * A2 + b[1] * MatrixType::Identity(A.rows(), A.cols()); @@ -98,7 +98,7 @@ template <typename MatrixType> void matrix_exp_pade7(const MatrixType &A, MatrixType &U, MatrixType &V) { typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar; - const RealScalar b[] = {17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.}; + const RealScalar b[] = {17297280.L, 8648640.L, 1995840.L, 277200.L, 25200.L, 1512.L, 56.L, 1.L}; const MatrixType A2 = A * A; const MatrixType A4 = A2 * A2; const MatrixType A6 = A4 * A2; @@ -118,8 +118,8 @@ template <typename MatrixType> void matrix_exp_pade9(const MatrixType &A, MatrixType &U, MatrixType &V) { typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar; - const RealScalar b[] = {17643225600., 8821612800., 2075673600., 302702400., 30270240., - 2162160., 110880., 3960., 90., 1.}; + const RealScalar b[] = {17643225600.L, 8821612800.L, 2075673600.L, 302702400.L, 30270240.L, + 2162160.L, 110880.L, 3960.L, 90.L, 1.L}; const MatrixType A2 = A * A; const MatrixType A4 = A2 * A2; const MatrixType A6 = A4 * A2; @@ -139,9 +139,9 @@ template <typename MatrixType> void matrix_exp_pade13(const MatrixType &A, MatrixType &U, MatrixType &V) { typedef typename NumTraits<typename traits<MatrixType>::Scalar>::Real RealScalar; - const RealScalar b[] = {64764752532480000., 32382376266240000., 7771770303897600., - 1187353796428800., 129060195264000., 10559470521600., 670442572800., - 33522128640., 1323241920., 40840800., 960960., 16380., 182., 1.}; + const RealScalar b[] = {64764752532480000.L, 32382376266240000.L, 7771770303897600.L, + 1187353796428800.L, 129060195264000.L, 10559470521600.L, 670442572800.L, + 33522128640.L, 1323241920.L, 40840800.L, 960960.L, 16380.L, 182.L, 1.L}; const MatrixType A2 = A * A; const MatrixType A4 = A2 * A2; const MatrixType A6 = A4 * A2; @@ -210,9 +210,9 @@ struct matrix_exp_computeUV<MatrixType, float> using std::pow; const float l1norm = arg.cwiseAbs().colwise().sum().maxCoeff(); squarings = 0; - if (l1norm < 4.258730016922831e-001) { + if (l1norm < 4.258730016922831e-001f) { matrix_exp_pade3(arg, U, V); - } else if (l1norm < 1.880152677804762e+000) { + } else if (l1norm < 1.880152677804762e+000f) { matrix_exp_pade5(arg, U, V); } else { const float maxnorm = 3.925724783138660f; diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h b/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h index 8f7a6f3b0..db2449d02 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h @@ -132,6 +132,7 @@ template <typename EivalsType, typename Cluster> void matrix_function_partition_eigenvalues(const EivalsType& eivals, std::list<Cluster>& clusters) { typedef typename EivalsType::Index Index; + typedef typename EivalsType::RealScalar RealScalar; for (Index i=0; i<eivals.rows(); ++i) { // Find cluster containing i-th ei'val, adding a new cluster if necessary typename std::list<Cluster>::iterator qi = matrix_function_find_cluster(i, clusters); @@ -145,7 +146,7 @@ void matrix_function_partition_eigenvalues(const EivalsType& eivals, std::list<C // Look for other element to add to the set for (Index j=i+1; j<eivals.rows(); ++j) { - if (abs(eivals(j) - eivals(i)) <= matrix_function_separation + if (abs(eivals(j) - eivals(i)) <= RealScalar(matrix_function_separation) && std::find(qi->begin(), qi->end(), j) == qi->end()) { typename std::list<Cluster>::iterator qj = matrix_function_find_cluster(j, clusters); if (qj == clusters.end()) { @@ -403,11 +404,10 @@ struct matrix_function_compute<MatrixType, 0> typedef internal::traits<MatrixType> Traits; typedef typename Traits::Scalar Scalar; static const int Rows = Traits::RowsAtCompileTime, Cols = Traits::ColsAtCompileTime; - static const int Options = MatrixType::Options; static const int MaxRows = Traits::MaxRowsAtCompileTime, MaxCols = Traits::MaxColsAtCompileTime; typedef std::complex<Scalar> ComplexScalar; - typedef Matrix<ComplexScalar, Rows, Cols, Options, MaxRows, MaxCols> ComplexMatrix; + typedef Matrix<ComplexScalar, Rows, Cols, 0, MaxRows, MaxCols> ComplexMatrix; ComplexMatrix CA = A.template cast<ComplexScalar>(); ComplexMatrix Cresult; @@ -508,9 +508,8 @@ template<typename Derived> class MatrixFunctionReturnValue typedef internal::traits<NestedEvalTypeClean> Traits; static const int RowsAtCompileTime = Traits::RowsAtCompileTime; static const int ColsAtCompileTime = Traits::ColsAtCompileTime; - static const int Options = NestedEvalTypeClean::Options; typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar; - typedef Matrix<ComplexScalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType; + typedef Matrix<ComplexScalar, Dynamic, Dynamic, 0, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType; typedef internal::MatrixFunctionAtomic<DynMatrixType> AtomicType; AtomicType atomic(m_f); diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h b/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h index e43e86e90..1acfbed9e 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h @@ -37,6 +37,7 @@ template <typename MatrixType> void matrix_log_compute_2x2(const MatrixType& A, MatrixType& result) { typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; using std::abs; using std::ceil; using std::imag; @@ -54,14 +55,14 @@ void matrix_log_compute_2x2(const MatrixType& A, MatrixType& result) { result(0,1) = A(0,1) / A(0,0); } - else if ((abs(A(0,0)) < 0.5*abs(A(1,1))) || (abs(A(0,0)) > 2*abs(A(1,1)))) + else if ((abs(A(0,0)) < RealScalar(0.5)*abs(A(1,1))) || (abs(A(0,0)) > 2*abs(A(1,1)))) { result(0,1) = A(0,1) * (logA11 - logA00) / y; } else { // computation in previous branch is inaccurate if A(1,1) \approx A(0,0) - int unwindingNumber = static_cast<int>(ceil((imag(logA11 - logA00) - EIGEN_PI) / (2*EIGEN_PI))); + int unwindingNumber = static_cast<int>(ceil((imag(logA11 - logA00) - RealScalar(EIGEN_PI)) / RealScalar(2*EIGEN_PI))); result(0,1) = A(0,1) * (numext::log1p(y/A(0,0)) + Scalar(0,2*EIGEN_PI*unwindingNumber)) / y; } } @@ -232,8 +233,8 @@ void matrix_log_compute_big(const MatrixType& A, MatrixType& result) MatrixType T = A, sqrtT; int maxPadeDegree = matrix_log_max_pade_degree<Scalar>::value; - const RealScalar maxNormForPade = maxPadeDegree<= 5? 5.3149729967117310e-1: // single precision - maxPadeDegree<= 7? 2.6429608311114350e-1: // double precision + const RealScalar maxNormForPade = maxPadeDegree<= 5? 5.3149729967117310e-1L: // single precision + maxPadeDegree<= 7? 2.6429608311114350e-1L: // double precision maxPadeDegree<= 8? 2.32777776523703892094e-1L: // extended precision maxPadeDegree<=10? 1.05026503471351080481093652651105e-1L: // double-double 1.1880960220216759245467951592883642e-1L; // quadruple precision @@ -333,9 +334,8 @@ public: typedef internal::traits<DerivedEvalTypeClean> Traits; static const int RowsAtCompileTime = Traits::RowsAtCompileTime; static const int ColsAtCompileTime = Traits::ColsAtCompileTime; - static const int Options = DerivedEvalTypeClean::Options; typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar; - typedef Matrix<ComplexScalar, Dynamic, Dynamic, Options, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType; + typedef Matrix<ComplexScalar, Dynamic, Dynamic, 0, RowsAtCompileTime, ColsAtCompileTime> DynMatrixType; typedef internal::MatrixLogarithmAtomic<DynMatrixType> AtomicType; AtomicType atomic; diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index f37d31c3f..ebc433d89 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -196,11 +196,11 @@ void MatrixPowerAtomic<MatrixType>::computeBig(ResultType& res) const { using std::ldexp; const int digits = std::numeric_limits<RealScalar>::digits; - const RealScalar maxNormForPade = digits <= 24? 4.3386528e-1f: // sigle precision - digits <= 53? 2.789358995219730e-1: // double precision - digits <= 64? 2.4471944416607995472e-1L: // extended precision - digits <= 106? 1.1016843812851143391275867258512e-1L: // double-double - 9.134603732914548552537150753385375e-2L; // quadruple precision + const RealScalar maxNormForPade = digits <= 24? 4.3386528e-1L // single precision + : digits <= 53? 2.789358995219730e-1L // double precision + : digits <= 64? 2.4471944416607995472e-1L // extended precision + : digits <= 106? 1.1016843812851143391275867258512e-1L // double-double + : 9.134603732914548552537150753385375e-2L; // quadruple precision MatrixType IminusT, sqrtT, T = m_A.template triangularView<Upper>(); RealScalar normIminusT; int degree, degree2, numberOfSquareRoots = 0; @@ -264,7 +264,7 @@ inline int MatrixPowerAtomic<MatrixType>::getPadeDegree(long double normIminusT) 1.999045567181744e-1L, 2.789358995219730e-1L }; #elif LDBL_MANT_DIG <= 64 const int maxPadeDegree = 8; - const double maxNormForPade[] = { 6.3854693117491799460e-3L /* degree = 3 */ , 2.6394893435456973676e-2L, + const long double maxNormForPade[] = { 6.3854693117491799460e-3L /* degree = 3 */ , 2.6394893435456973676e-2L, 6.4216043030404063729e-2L, 1.1701165502926694307e-1L, 1.7904284231268670284e-1L, 2.4471944416607995472e-1L }; #elif LDBL_MANT_DIG <= 106 const int maxPadeDegree = 10; @@ -298,7 +298,7 @@ MatrixPowerAtomic<MatrixType>::computeSuperDiag(const ComplexScalar& curr, const ComplexScalar logCurr = log(curr); ComplexScalar logPrev = log(prev); - int unwindingNumber = ceil((numext::imag(logCurr - logPrev) - EIGEN_PI) / (2*EIGEN_PI)); + int unwindingNumber = ceil((numext::imag(logCurr - logPrev) - RealScalar(EIGEN_PI)) / RealScalar(2*EIGEN_PI)); ComplexScalar w = numext::log1p((curr-prev)/prev)/RealScalar(2) + ComplexScalar(0, EIGEN_PI*unwindingNumber); return RealScalar(2) * exp(RealScalar(0.5) * p * (logCurr + logPrev)) * sinh(p * w) / (curr - prev); } @@ -383,7 +383,7 @@ class MatrixPower : internal::noncopyable private: typedef std::complex<RealScalar> ComplexScalar; - typedef Matrix<ComplexScalar, Dynamic, Dynamic, MatrixType::Options, + typedef Matrix<ComplexScalar, Dynamic, Dynamic, 0, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime> ComplexMatrix; /** \brief Reference to the base of matrix power. */ diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h b/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h index 9f08c6162..afd88ec4d 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h @@ -65,21 +65,6 @@ void matrix_sqrt_quasi_triangular_2x1_off_diagonal_block(const MatrixType& T, ty sqrtT.template block<2,1>(i,j) = A.fullPivLu().solve(rhs); } -// similar to compute1x1offDiagonalBlock() -template <typename MatrixType, typename ResultType> -void matrix_sqrt_quasi_triangular_2x2_off_diagonal_block(const MatrixType& T, typename MatrixType::Index i, typename MatrixType::Index j, ResultType& sqrtT) -{ - typedef typename traits<MatrixType>::Scalar Scalar; - Matrix<Scalar,2,2> A = sqrtT.template block<2,2>(i,i); - Matrix<Scalar,2,2> B = sqrtT.template block<2,2>(j,j); - Matrix<Scalar,2,2> C = T.template block<2,2>(i,j); - if (j-i > 2) - C -= sqrtT.block(i, i+2, 2, j-i-2) * sqrtT.block(i+2, j, j-i-2, 2); - Matrix<Scalar,2,2> X; - matrix_sqrt_quasi_triangular_solve_auxiliary_equation(X, A, B, C); - sqrtT.template block<2,2>(i,j) = X; -} - // solves the equation A X + X B = C where all matrices are 2-by-2 template <typename MatrixType> void matrix_sqrt_quasi_triangular_solve_auxiliary_equation(MatrixType& X, const MatrixType& A, const MatrixType& B, const MatrixType& C) @@ -98,13 +83,13 @@ void matrix_sqrt_quasi_triangular_solve_auxiliary_equation(MatrixType& X, const coeffMatrix.coeffRef(2,3) = B.coeff(1,0); coeffMatrix.coeffRef(3,1) = A.coeff(1,0); coeffMatrix.coeffRef(3,2) = B.coeff(0,1); - + Matrix<Scalar,4,1> rhs; rhs.coeffRef(0) = C.coeff(0,0); rhs.coeffRef(1) = C.coeff(0,1); rhs.coeffRef(2) = C.coeff(1,0); rhs.coeffRef(3) = C.coeff(1,1); - + Matrix<Scalar,4,1> result; result = coeffMatrix.fullPivLu().solve(rhs); @@ -114,6 +99,20 @@ void matrix_sqrt_quasi_triangular_solve_auxiliary_equation(MatrixType& X, const X.coeffRef(1,1) = result.coeff(3); } +// similar to compute1x1offDiagonalBlock() +template <typename MatrixType, typename ResultType> +void matrix_sqrt_quasi_triangular_2x2_off_diagonal_block(const MatrixType& T, typename MatrixType::Index i, typename MatrixType::Index j, ResultType& sqrtT) +{ + typedef typename traits<MatrixType>::Scalar Scalar; + Matrix<Scalar,2,2> A = sqrtT.template block<2,2>(i,i); + Matrix<Scalar,2,2> B = sqrtT.template block<2,2>(j,j); + Matrix<Scalar,2,2> C = T.template block<2,2>(i,j); + if (j-i > 2) + C -= sqrtT.block(i, i+2, 2, j-i-2) * sqrtT.block(i+2, j, j-i-2, 2); + Matrix<Scalar,2,2> X; + matrix_sqrt_quasi_triangular_solve_auxiliary_equation(X, A, B, C); + sqrtT.template block<2,2>(i,j) = X; +} // pre: T is quasi-upper-triangular and sqrtT is a zero matrix of the same size // post: the diagonal blocks of sqrtT are the square roots of the diagonal blocks of T |