1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_JACOBI_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_JACOBI_H 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup Jacobi_Module 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \jacobi_module 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class JacobiRotation 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Rotation given by a cosine-sine pair. 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class represents a Jacobi or Givens rotation. 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This is a 2D rotation in the plane \c J of angle \f$ \theta \f$ defined by 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * its cosine \c c and sine \c s as follow: 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ J = \left ( \begin{array}{cc} c & \overline s \\ -s & \overline c \end{array} \right ) \f$ 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * You can apply the respective counter-clockwise rotation to a column vector \c v by 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * applying its adjoint on the left: \f$ v = J^* v \f$ that translates to the following Eigen code: 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \code 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * v.applyOnTheLeft(J.adjoint()); 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \endcode 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> class JacobiRotation 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<Scalar>::Real RealScalar; 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Default constructor without any initialization. */ 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiRotation() {} 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Construct a planar rotation from a cosine-sine pair (\a c, \c s). */ 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiRotation(const Scalar& c, const Scalar& s) : m_c(c), m_s(s) {} 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar& c() { return m_c; } 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar c() const { return m_c; } 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar& s() { return m_s; } 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar s() const { return m_s; } 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Concatenates two planar rotation */ 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiRotation operator*(const JacobiRotation& other) 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using numext::conj; 547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return JacobiRotation(m_c * other.m_c - conj(m_s) * other.m_s, 557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez conj(m_c * conj(other.m_s) + conj(m_s) * conj(other.m_c))); 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Returns the transposed transformation */ 597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez JacobiRotation transpose() const { using numext::conj; return JacobiRotation(m_c, -conj(m_s)); } 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Returns the adjoint transformation */ 627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez JacobiRotation adjoint() const { using numext::conj; return JacobiRotation(conj(m_c), -m_s); } 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Derived> 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool makeJacobi(const MatrixBase<Derived>&, typename Derived::Index p, typename Derived::Index q); 667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool makeJacobi(const RealScalar& x, const Scalar& y, const RealScalar& z); 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void makeGivens(const Scalar& p, const Scalar& q, Scalar* z=0); 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void makeGivens(const Scalar& p, const Scalar& q, Scalar* z, internal::true_type); 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void makeGivens(const Scalar& p, const Scalar& q, Scalar* z, internal::false_type); 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar m_c, m_s; 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** Makes \c *this as a Jacobi rotation \a J such that applying \a J on both the right and left sides of the selfadjoint 2x2 matrix 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ B = \left ( \begin{array}{cc} x & y \\ \overline y & z \end{array} \right )\f$ yields a diagonal matrix \f$ A = J^* B J \f$ 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixBase::makeJacobi(const MatrixBase<Derived>&, Index, Index), MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> 837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezbool JacobiRotation<Scalar>::makeJacobi(const RealScalar& x, const Scalar& y, const RealScalar& z) 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::sqrt; 867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::abs; 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<Scalar>::Real RealScalar; 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(y == Scalar(0)) 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_c = Scalar(1); 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_s = Scalar(0); 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return false; 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar tau = (x-z)/(RealScalar(2)*abs(y)); 977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar w = sqrt(numext::abs2(tau) + RealScalar(1)); 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar t; 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(tau>RealScalar(0)) 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath t = RealScalar(1) / (tau + w); 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath t = RealScalar(1) / (tau - w); 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar sign_t = t > RealScalar(0) ? RealScalar(1) : RealScalar(-1); 1087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar n = RealScalar(1) / sqrt(numext::abs2(t)+RealScalar(1)); 1097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_s = - sign_t * (numext::conj(y) / abs(y)) * abs(t) * n; 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_c = n; 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return true; 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** Makes \c *this as a Jacobi rotation \c J such that applying \a J on both the right and left sides of the 2x2 selfadjoint matrix 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ B = \left ( \begin{array}{cc} \text{this}_{pp} & \text{this}_{pq} \\ (\text{this}_{pq})^* & \text{this}_{qq} \end{array} \right )\f$ yields 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * a diagonal matrix \f$ A = J^* B J \f$ 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include Jacobi_makeJacobi.cpp 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude Jacobi_makeJacobi.out 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa JacobiRotation::makeJacobi(RealScalar, Scalar, RealScalar), MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathinline bool JacobiRotation<Scalar>::makeJacobi(const MatrixBase<Derived>& m, typename Derived::Index p, typename Derived::Index q) 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return makeJacobi(numext::real(m.coeff(p,p)), m.coeff(p,q), numext::real(m.coeff(q,q))); 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** Makes \c *this as a Givens rotation \c G such that applying \f$ G^* \f$ to the left of the vector 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ V = \left ( \begin{array}{c} p \\ q \end{array} \right )\f$ yields: 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ G^* V = \left ( \begin{array}{c} r \\ 0 \end{array} \right )\f$. 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The value of \a z is returned if \a z is not null (the default is null). 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Also note that G is built such that the cosine is always real. 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include Jacobi_makeGivens.cpp 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude Jacobi_makeGivens.out 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This function implements the continuous Givens rotation generation algorithm 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * found in Anderson (2000), Discontinuous Plane Rotations and the Symmetric Eigenvalue Problem. 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * LAPACK Working Note 150, University of Tennessee, UT-CS-00-454, December 4, 2000. 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid JacobiRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* z) 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath makeGivens(p, q, z, typename internal::conditional<NumTraits<Scalar>::IsComplex, internal::true_type, internal::false_type>::type()); 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// specialization for complexes 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid JacobiRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* r, internal::true_type) 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::sqrt; 1597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::abs; 1607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using numext::conj; 1617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(q==Scalar(0)) 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_c = numext::real(p)<0 ? Scalar(-1) : Scalar(1); 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_s = 0; 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(r) *r = m_c * p; 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(p==Scalar(0)) 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_c = 0; 1717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_s = -q/abs(q); 1727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(r) *r = abs(q); 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar p1 = numext::norm1(p); 1777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar q1 = numext::norm1(q); 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(p1>=q1) 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar ps = p / p1; 1817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar p2 = numext::abs2(ps); 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar qs = q / p1; 1837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar q2 = numext::abs2(qs); 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar u = sqrt(RealScalar(1) + q2/p2); 1867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(numext::real(p)<RealScalar(0)) 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath u = -u; 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_c = Scalar(1)/u; 1907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_s = -qs*conj(ps)*(m_c/p2); 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(r) *r = p * u; 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar ps = p / q1; 1967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar p2 = numext::abs2(ps); 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar qs = q / q1; 1987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar q2 = numext::abs2(qs); 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar u = q1 * sqrt(p2 + q2); 2017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(numext::real(p)<RealScalar(0)) 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath u = -u; 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez p1 = abs(p); 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ps = p/p1; 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_c = p1/u; 2077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_s = -conj(ps) * (q/u); 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(r) *r = ps * u; 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// specialization for reals 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid JacobiRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* r, internal::false_type) 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 2177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::sqrt; 2187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::abs; 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(q==Scalar(0)) 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_c = p<Scalar(0) ? Scalar(-1) : Scalar(1); 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_s = Scalar(0); 2237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(r) *r = abs(p); 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(p==Scalar(0)) 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_c = Scalar(0); 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_s = q<Scalar(0) ? Scalar(1) : Scalar(-1); 2297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(r) *r = abs(q); 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 2317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if(abs(p) > abs(q)) 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar t = q/p; 2347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar u = sqrt(Scalar(1) + numext::abs2(t)); 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(p<Scalar(0)) 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath u = -u; 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_c = Scalar(1)/u; 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_s = -t * m_c; 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(r) *r = p * u; 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar t = p/q; 2447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar u = sqrt(Scalar(1) + numext::abs2(t)); 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(q<Scalar(0)) 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath u = -u; 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_s = -Scalar(1)/u; 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_c = -t * m_s; 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(r) *r = q * u; 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/**************************************************************************************** 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath* Implementation of MatrixBase methods 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath****************************************************************************************/ 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \jacobi_module 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Applies the clock wise 2D rotation \a j to the set of 2D vectors of cordinates \a x and \a y: 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ \left ( \begin{array}{cc} x \\ y \end{array} \right ) = J \left ( \begin{array}{cc} x \\ y \end{array} \right ) \f$ 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename VectorX, typename VectorY, typename OtherScalar> 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, const JacobiRotation<OtherScalar>& j); 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \jacobi_module 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Applies the rotation in the plane \a j to the rows \a p and \a q of \c *this, i.e., it computes B = J * B, 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * with \f$ B = \left ( \begin{array}{cc} \text{*this.row}(p) \\ \text{*this.row}(q) \end{array} \right ) \f$. 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa class JacobiRotation, MatrixBase::applyOnTheRight(), internal::apply_rotation_in_the_plane() 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename OtherScalar> 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathinline void MatrixBase<Derived>::applyOnTheLeft(Index p, Index q, const JacobiRotation<OtherScalar>& j) 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RowXpr x(this->row(p)); 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RowXpr y(this->row(q)); 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::apply_rotation_in_the_plane(x, y, j); 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup Jacobi_Module 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Applies the rotation in the plane \a j to the columns \a p and \a q of \c *this, i.e., it computes B = B * J 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * with \f$ B = \left ( \begin{array}{cc} \text{*this.col}(p) & \text{*this.col}(q) \end{array} \right ) \f$. 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa class JacobiRotation, MatrixBase::applyOnTheLeft(), internal::apply_rotation_in_the_plane() 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename OtherScalar> 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathinline void MatrixBase<Derived>::applyOnTheRight(Index p, Index q, const JacobiRotation<OtherScalar>& j) 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColXpr x(this->col(p)); 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColXpr y(this->col(q)); 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::apply_rotation_in_the_plane(x, y, j.transpose()); 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename VectorX, typename VectorY, typename OtherScalar> 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid /*EIGEN_DONT_INLINE*/ apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, const JacobiRotation<OtherScalar>& j) 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename VectorX::Index Index; 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename VectorX::Scalar Scalar; 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { PacketSize = packet_traits<Scalar>::size }; 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename packet_traits<Scalar>::type Packet; 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(_x.size() == _y.size()); 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index size = _x.size(); 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index incrx = _x.innerStride(); 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index incry = _y.innerStride(); 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* EIGEN_RESTRICT x = &_x.coeffRef(0); 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* EIGEN_RESTRICT y = &_y.coeffRef(0); 3147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez OtherScalar c = j.c(); 3167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez OtherScalar s = j.s(); 3177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (c==OtherScalar(1) && s==OtherScalar(0)) 3187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return; 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /*** dynamic-size vectorized paths ***/ 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(VectorX::SizeAtCompileTime == Dynamic && 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath (VectorX::Flags & VectorY::Flags & PacketAccessBit) && 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ((incrx==1 && incry==1) || PacketSize == 1)) 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // both vectors are sequentially stored in memory => vectorization 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { Peeling = 2 }; 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index alignedStart = internal::first_aligned(y, size); 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index alignedEnd = alignedStart + ((size-alignedStart)/PacketSize)*PacketSize; 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const Packet pc = pset1<Packet>(c); 3337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const Packet ps = pset1<Packet>(s); 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath conj_helper<Packet,Packet,NumTraits<Scalar>::IsComplex,false> pcj; 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index i=0; i<alignedStart; ++i) 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar xi = x[i]; 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar yi = y[i]; 3407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x[i] = c * xi + numext::conj(s) * yi; 3417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez y[i] = -s * xi + numext::conj(c) * yi; 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* EIGEN_RESTRICT px = x + alignedStart; 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* EIGEN_RESTRICT py = y + alignedStart; 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(internal::first_aligned(x, size)==alignedStart) 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index i=alignedStart; i<alignedEnd; i+=PacketSize) 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Packet xi = pload<Packet>(px); 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Packet yi = pload<Packet>(py); 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pstore(px, padd(pmul(pc,xi),pcj.pmul(ps,yi))); 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pstore(py, psub(pcj.pmul(pc,yi),pmul(ps,xi))); 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath px += PacketSize; 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath py += PacketSize; 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index peelingEnd = alignedStart + ((size-alignedStart)/(Peeling*PacketSize))*(Peeling*PacketSize); 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index i=alignedStart; i<peelingEnd; i+=Peeling*PacketSize) 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Packet xi = ploadu<Packet>(px); 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Packet xi1 = ploadu<Packet>(px+PacketSize); 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Packet yi = pload <Packet>(py); 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Packet yi1 = pload <Packet>(py+PacketSize); 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pstoreu(px, padd(pmul(pc,xi),pcj.pmul(ps,yi))); 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pstoreu(px+PacketSize, padd(pmul(pc,xi1),pcj.pmul(ps,yi1))); 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pstore (py, psub(pcj.pmul(pc,yi),pmul(ps,xi))); 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pstore (py+PacketSize, psub(pcj.pmul(pc,yi1),pmul(ps,xi1))); 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath px += Peeling*PacketSize; 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath py += Peeling*PacketSize; 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(alignedEnd!=peelingEnd) 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Packet xi = ploadu<Packet>(x+peelingEnd); 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Packet yi = pload <Packet>(y+peelingEnd); 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pstoreu(x+peelingEnd, padd(pmul(pc,xi),pcj.pmul(ps,yi))); 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pstore (y+peelingEnd, psub(pcj.pmul(pc,yi),pmul(ps,xi))); 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index i=alignedEnd; i<size; ++i) 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar xi = x[i]; 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar yi = y[i]; 3887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x[i] = c * xi + numext::conj(s) * yi; 3897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez y[i] = -s * xi + numext::conj(c) * yi; 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /*** fixed-size vectorized path ***/ 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(VectorX::SizeAtCompileTime != Dynamic && 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath (VectorX::Flags & VectorY::Flags & PacketAccessBit) && 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath (VectorX::Flags & VectorY::Flags & AlignedBit)) 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 3987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const Packet pc = pset1<Packet>(c); 3997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const Packet ps = pset1<Packet>(s); 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath conj_helper<Packet,Packet,NumTraits<Scalar>::IsComplex,false> pcj; 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* EIGEN_RESTRICT px = x; 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* EIGEN_RESTRICT py = y; 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index i=0; i<size; i+=PacketSize) 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Packet xi = pload<Packet>(px); 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Packet yi = pload<Packet>(py); 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pstore(px, padd(pmul(pc,xi),pcj.pmul(ps,yi))); 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pstore(py, psub(pcj.pmul(pc,yi),pmul(ps,xi))); 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath px += PacketSize; 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath py += PacketSize; 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /*** non-vectorized path ***/ 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index i=0; i<size; ++i) 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar xi = *x; 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar yi = *y; 4217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *x = c * xi + numext::conj(s) * yi; 4227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *y = -s * xi + numext::conj(c) * yi; 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x += incrx; 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath y += incry; 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_JACOBI_H 434