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 (at) gmail.com> 5 // Copyright (C) 2013-2016 Gael Guennebaud <gael.guennebaud (at) 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 56
