![]() |
Eigen
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2009-2010 Benoit Jacob <jacob.benoit.1@gmail.com> 00005 // Copyright (C) 2013-2016 Gael Guennebaud <gael.guennebaud@inria.fr> 00006 // 00007 // This Source Code Form is subject to the terms of the Mozilla 00008 // Public License v. 2.0. If a copy of the MPL was not distributed 00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00010 00011 #ifndef EIGEN_REALSVD2X2_H 00012 #define EIGEN_REALSVD2X2_H 00013 00014 namespace Eigen { 00015 00016 namespace internal { 00017 00018 template<typename MatrixType, typename RealScalar, typename Index> 00019 void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, 00020 JacobiRotation<RealScalar> *j_left, 00021 JacobiRotation<RealScalar> *j_right) 00022 { 00023 using std::sqrt; 00024 using std::abs; 00025 Matrix<RealScalar,2,2> m; 00026 m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)), 00027 numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q)); 00028 JacobiRotation<RealScalar> rot1; 00029 RealScalar t = m.coeff(0,0) + m.coeff(1,1); 00030 RealScalar d = m.coeff(1,0) - m.coeff(0,1); 00031 00032 if(abs(d) < (std::numeric_limits<RealScalar>::min)()) 00033 { 00034 rot1.s() = RealScalar(0); 00035 rot1.c() = RealScalar(1); 00036 } 00037 else 00038 { 00039 // If d!=0, then t/d cannot overflow because the magnitude of the 00040 // entries forming d are not too small compared to the ones forming t. 00041 RealScalar u = t / d; 00042 RealScalar tmp = sqrt(RealScalar(1) + numext::abs2(u)); 00043 rot1.s() = RealScalar(1) / tmp; 00044 rot1.c() = u / tmp; 00045 } 00046 m.applyOnTheLeft(0,1,rot1); 00047 j_right->makeJacobi(m,0,1); 00048 *j_left = rot1 * j_right->transpose(); 00049 } 00050 00051 } // end namespace internal 00052 00053 } // end namespace Eigen 00054 00055 #endif // EIGEN_REALSVD2X2_H