| Austin Schuh | 189376f | 2018-12-20 22:11:15 +1100 | [diff] [blame^] | 1 | // This file is part of Eigen, a lightweight C++ template library |
| 2 | // for linear algebra. |
| 3 | // |
| 4 | // Copyright (C) 2009-2010 Benoit Jacob <jacob.benoit.1@gmail.com> |
| 5 | // Copyright (C) 2013-2016 Gael Guennebaud <gael.guennebaud@inria.fr> |
| 6 | // |
| 7 | // This Source Code Form is subject to the terms of the Mozilla |
| 8 | // Public License v. 2.0. If a copy of the MPL was not distributed |
| 9 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
| 10 | |
| 11 | #ifndef EIGEN_REALSVD2X2_H |
| 12 | #define EIGEN_REALSVD2X2_H |
| 13 | |
| 14 | namespace Eigen { |
| 15 | |
| 16 | namespace internal { |
| 17 | |
| 18 | template<typename MatrixType, typename RealScalar, typename Index> |
| 19 | void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, |
| 20 | JacobiRotation<RealScalar> *j_left, |
| 21 | JacobiRotation<RealScalar> *j_right) |
| 22 | { |
| 23 | using std::sqrt; |
| 24 | using std::abs; |
| 25 | Matrix<RealScalar,2,2> m; |
| 26 | m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)), |
| 27 | numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q)); |
| 28 | JacobiRotation<RealScalar> rot1; |
| 29 | RealScalar t = m.coeff(0,0) + m.coeff(1,1); |
| 30 | RealScalar d = m.coeff(1,0) - m.coeff(0,1); |
| 31 | |
| 32 | if(abs(d) < (std::numeric_limits<RealScalar>::min)()) |
| 33 | { |
| 34 | rot1.s() = RealScalar(0); |
| 35 | rot1.c() = RealScalar(1); |
| 36 | } |
| 37 | else |
| 38 | { |
| 39 | // If d!=0, then t/d cannot overflow because the magnitude of the |
| 40 | // entries forming d are not too small compared to the ones forming t. |
| 41 | RealScalar u = t / d; |
| 42 | RealScalar tmp = sqrt(RealScalar(1) + numext::abs2(u)); |
| 43 | rot1.s() = RealScalar(1) / tmp; |
| 44 | rot1.c() = u / tmp; |
| 45 | } |
| 46 | m.applyOnTheLeft(0,1,rot1); |
| 47 | j_right->makeJacobi(m,0,1); |
| 48 | *j_left = rot1 * j_right->transpose(); |
| 49 | } |
| 50 | |
| 51 | } // end namespace internal |
| 52 | |
| 53 | } // end namespace Eigen |
| 54 | |
| 55 | #endif // EIGEN_REALSVD2X2_H |