12b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// This file is part of Eigen, a lightweight C++ template library 22b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// for linear algebra. 32b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// 42b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Copyright (C) 2009-2010 Benoit Jacob <jacob.benoit.1@gmail.com> 52b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Copyright (C) 2013-2016 Gael Guennebaud <gael.guennebaud@inria.fr> 62b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// 72b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// This Source Code Form is subject to the terms of the Mozilla 82b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Public License v. 2.0. If a copy of the MPL was not distributed 92b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#ifndef EIGEN_REALSVD2X2_H 122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#define EIGEN_REALSVD2X2_H 132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangnamespace Eigen { 152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangnamespace internal { 172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename MatrixType, typename RealScalar, typename Index> 192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangvoid real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, 202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang JacobiRotation<RealScalar> *j_left, 212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang JacobiRotation<RealScalar> *j_right) 222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang{ 232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang using std::sqrt; 242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang using std::abs; 252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Matrix<RealScalar,2,2> m; 262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)), 272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q)); 282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang JacobiRotation<RealScalar> rot1; 292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang RealScalar t = m.coeff(0,0) + m.coeff(1,1); 302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang RealScalar d = m.coeff(1,0) - m.coeff(0,1); 312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(abs(d) < (std::numeric_limits<RealScalar>::min)()) 332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang rot1.s() = RealScalar(0); 352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang rot1.c() = RealScalar(1); 362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang else 382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // If d!=0, then t/d cannot overflow because the magnitude of the 402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // entries forming d are not too small compared to the ones forming t. 412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang RealScalar u = t / d; 422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang RealScalar tmp = sqrt(RealScalar(1) + numext::abs2(u)); 432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang rot1.s() = RealScalar(1) / tmp; 442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang rot1.c() = u / tmp; 452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m.applyOnTheLeft(0,1,rot1); 472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang j_right->makeJacobi(m,0,1); 482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang *j_left = rot1 * j_right->transpose(); 492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} 502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} // end namespace internal 522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang} // end namespace Eigen 542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang#endif // EIGEN_REALSVD2X2_H 56