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 14namespace Eigen { 15 16namespace internal { 17 18template<typename MatrixType, typename RealScalar, typename Index> 19void 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