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 { 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return JacobiRotation(m_c * other.m_c - internal::conj(m_s) * other.m_s, 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::conj(m_c * internal::conj(other.m_s) + internal::conj(m_s) * internal::conj(other.m_c))); 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Returns the transposed transformation */ 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiRotation transpose() const { return JacobiRotation(m_c, -internal::conj(m_s)); } 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Returns the adjoint transformation */ 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiRotation adjoint() const { return JacobiRotation(internal::conj(m_c), -m_s); } 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Derived> 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool makeJacobi(const MatrixBase<Derived>&, typename Derived::Index p, typename Derived::Index q); 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool makeJacobi(RealScalar x, Scalar y, RealScalar z); 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void makeGivens(const Scalar& p, const Scalar& q, Scalar* z=0); 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void makeGivens(const Scalar& p, const Scalar& q, Scalar* z, internal::true_type); 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void makeGivens(const Scalar& p, const Scalar& q, Scalar* z, internal::false_type); 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar m_c, m_s; 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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$ 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixBase::makeJacobi(const MatrixBase<Derived>&, Index, Index), MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathbool JacobiRotation<Scalar>::makeJacobi(RealScalar x, Scalar y, RealScalar z) 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<Scalar>::Real RealScalar; 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(y == Scalar(0)) 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_c = Scalar(1); 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_s = Scalar(0); 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return false; 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar tau = (x-z)/(RealScalar(2)*internal::abs(y)); 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar w = internal::sqrt(internal::abs2(tau) + RealScalar(1)); 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar t; 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(tau>RealScalar(0)) 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath t = RealScalar(1) / (tau + w); 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath t = RealScalar(1) / (tau - w); 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar sign_t = t > RealScalar(0) ? RealScalar(1) : RealScalar(-1); 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar n = RealScalar(1) / internal::sqrt(internal::abs2(t)+RealScalar(1)); 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_s = - sign_t * (internal::conj(y) / internal::abs(y)) * internal::abs(t) * n; 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_c = n; 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return true; 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ B = \left ( \begin{array}{cc} \text{this}_{pp} & \text{this}_{pq} \\ (\text{this}_{pq})^* & \text{this}_{qq} \end{array} \right )\f$ yields 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * a diagonal matrix \f$ A = J^* B J \f$ 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include Jacobi_makeJacobi.cpp 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude Jacobi_makeJacobi.out 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa JacobiRotation::makeJacobi(RealScalar, Scalar, RealScalar), MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathinline bool JacobiRotation<Scalar>::makeJacobi(const MatrixBase<Derived>& m, typename Derived::Index p, typename Derived::Index q) 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return makeJacobi(internal::real(m.coeff(p,p)), m.coeff(p,q), internal::real(m.coeff(q,q))); 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** Makes \c *this as a Givens rotation \c G such that applying \f$ G^* \f$ to the left of the vector 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ V = \left ( \begin{array}{c} p \\ q \end{array} \right )\f$ yields: 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ G^* V = \left ( \begin{array}{c} r \\ 0 \end{array} \right )\f$. 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The value of \a z is returned if \a z is not null (the default is null). 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Also note that G is built such that the cosine is always real. 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include Jacobi_makeGivens.cpp 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude Jacobi_makeGivens.out 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This function implements the continuous Givens rotation generation algorithm 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * found in Anderson (2000), Discontinuous Plane Rotations and the Symmetric Eigenvalue Problem. 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * LAPACK Working Note 150, University of Tennessee, UT-CS-00-454, December 4, 2000. 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid JacobiRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* z) 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath makeGivens(p, q, z, typename internal::conditional<NumTraits<Scalar>::IsComplex, internal::true_type, internal::false_type>::type()); 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// specialization for complexes 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid JacobiRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* r, internal::true_type) 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(q==Scalar(0)) 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_c = internal::real(p)<0 ? Scalar(-1) : Scalar(1); 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_s = 0; 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(r) *r = m_c * p; 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(p==Scalar(0)) 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_c = 0; 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_s = -q/internal::abs(q); 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(r) *r = internal::abs(q); 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar p1 = internal::norm1(p); 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar q1 = internal::norm1(q); 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(p1>=q1) 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar ps = p / p1; 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar p2 = internal::abs2(ps); 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar qs = q / p1; 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar q2 = internal::abs2(qs); 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar u = internal::sqrt(RealScalar(1) + q2/p2); 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(internal::real(p)<RealScalar(0)) 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath u = -u; 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_c = Scalar(1)/u; 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_s = -qs*internal::conj(ps)*(m_c/p2); 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(r) *r = p * u; 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar ps = p / q1; 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar p2 = internal::abs2(ps); 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar qs = q / q1; 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar q2 = internal::abs2(qs); 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar u = q1 * internal::sqrt(p2 + q2); 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(internal::real(p)<RealScalar(0)) 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath u = -u; 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath p1 = internal::abs(p); 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ps = p/p1; 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_c = p1/u; 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_s = -internal::conj(ps) * (q/u); 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(r) *r = ps * u; 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// specialization for reals 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid JacobiRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* r, internal::false_type) 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(q==Scalar(0)) 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_c = p<Scalar(0) ? Scalar(-1) : Scalar(1); 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_s = Scalar(0); 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(r) *r = internal::abs(p); 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(p==Scalar(0)) 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_c = Scalar(0); 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_s = q<Scalar(0) ? Scalar(1) : Scalar(-1); 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(r) *r = internal::abs(q); 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(internal::abs(p) > internal::abs(q)) 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar t = q/p; 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar u = internal::sqrt(Scalar(1) + internal::abs2(t)); 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(p<Scalar(0)) 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath u = -u; 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_c = Scalar(1)/u; 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_s = -t * m_c; 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(r) *r = p * u; 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar t = p/q; 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar u = internal::sqrt(Scalar(1) + internal::abs2(t)); 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(q<Scalar(0)) 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath u = -u; 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_s = -Scalar(1)/u; 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_c = -t * m_s; 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(r) *r = q * u; 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/**************************************************************************************** 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath* Implementation of MatrixBase methods 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath****************************************************************************************/ 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \jacobi_module 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Applies the clock wise 2D rotation \a j to the set of 2D vectors of cordinates \a x and \a y: 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ \left ( \begin{array}{cc} x \\ y \end{array} \right ) = J \left ( \begin{array}{cc} x \\ y \end{array} \right ) \f$ 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename VectorX, typename VectorY, typename OtherScalar> 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, const JacobiRotation<OtherScalar>& j); 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \jacobi_module 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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, 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * with \f$ B = \left ( \begin{array}{cc} \text{*this.row}(p) \\ \text{*this.row}(q) \end{array} \right ) \f$. 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa class JacobiRotation, MatrixBase::applyOnTheRight(), internal::apply_rotation_in_the_plane() 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename OtherScalar> 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathinline void MatrixBase<Derived>::applyOnTheLeft(Index p, Index q, const JacobiRotation<OtherScalar>& j) 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RowXpr x(this->row(p)); 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RowXpr y(this->row(q)); 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::apply_rotation_in_the_plane(x, y, j); 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup Jacobi_Module 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * with \f$ B = \left ( \begin{array}{cc} \text{*this.col}(p) & \text{*this.col}(q) \end{array} \right ) \f$. 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa class JacobiRotation, MatrixBase::applyOnTheLeft(), internal::apply_rotation_in_the_plane() 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename OtherScalar> 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathinline void MatrixBase<Derived>::applyOnTheRight(Index p, Index q, const JacobiRotation<OtherScalar>& j) 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColXpr x(this->col(p)); 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColXpr y(this->col(q)); 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::apply_rotation_in_the_plane(x, y, j.transpose()); 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename VectorX, typename VectorY, typename OtherScalar> 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid /*EIGEN_DONT_INLINE*/ apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, const JacobiRotation<OtherScalar>& j) 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename VectorX::Index Index; 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename VectorX::Scalar Scalar; 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { PacketSize = packet_traits<Scalar>::size }; 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename packet_traits<Scalar>::type Packet; 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(_x.size() == _y.size()); 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index size = _x.size(); 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index incrx = _x.innerStride(); 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index incry = _y.innerStride(); 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* EIGEN_RESTRICT x = &_x.coeffRef(0); 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* EIGEN_RESTRICT y = &_y.coeffRef(0); 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /*** dynamic-size vectorized paths ***/ 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(VectorX::SizeAtCompileTime == Dynamic && 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath (VectorX::Flags & VectorY::Flags & PacketAccessBit) && 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ((incrx==1 && incry==1) || PacketSize == 1)) 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // both vectors are sequentially stored in memory => vectorization 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { Peeling = 2 }; 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index alignedStart = internal::first_aligned(y, size); 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index alignedEnd = alignedStart + ((size-alignedStart)/PacketSize)*PacketSize; 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Packet pc = pset1<Packet>(j.c()); 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Packet ps = pset1<Packet>(j.s()); 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath conj_helper<Packet,Packet,NumTraits<Scalar>::IsComplex,false> pcj; 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index i=0; i<alignedStart; ++i) 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar xi = x[i]; 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar yi = y[i]; 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x[i] = j.c() * xi + conj(j.s()) * yi; 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath y[i] = -j.s() * xi + conj(j.c()) * yi; 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* EIGEN_RESTRICT px = x + alignedStart; 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* EIGEN_RESTRICT py = y + alignedStart; 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(internal::first_aligned(x, size)==alignedStart) 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index i=alignedStart; i<alignedEnd; i+=PacketSize) 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Packet xi = pload<Packet>(px); 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Packet yi = pload<Packet>(py); 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pstore(px, padd(pmul(pc,xi),pcj.pmul(ps,yi))); 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pstore(py, psub(pcj.pmul(pc,yi),pmul(ps,xi))); 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath px += PacketSize; 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath py += PacketSize; 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index peelingEnd = alignedStart + ((size-alignedStart)/(Peeling*PacketSize))*(Peeling*PacketSize); 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index i=alignedStart; i<peelingEnd; i+=Peeling*PacketSize) 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Packet xi = ploadu<Packet>(px); 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Packet xi1 = ploadu<Packet>(px+PacketSize); 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Packet yi = pload <Packet>(py); 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Packet yi1 = pload <Packet>(py+PacketSize); 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pstoreu(px, padd(pmul(pc,xi),pcj.pmul(ps,yi))); 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pstoreu(px+PacketSize, padd(pmul(pc,xi1),pcj.pmul(ps,yi1))); 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pstore (py, psub(pcj.pmul(pc,yi),pmul(ps,xi))); 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pstore (py+PacketSize, psub(pcj.pmul(pc,yi1),pmul(ps,xi1))); 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath px += Peeling*PacketSize; 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath py += Peeling*PacketSize; 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(alignedEnd!=peelingEnd) 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Packet xi = ploadu<Packet>(x+peelingEnd); 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Packet yi = pload <Packet>(y+peelingEnd); 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pstoreu(x+peelingEnd, padd(pmul(pc,xi),pcj.pmul(ps,yi))); 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pstore (y+peelingEnd, psub(pcj.pmul(pc,yi),pmul(ps,xi))); 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index i=alignedEnd; i<size; ++i) 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar xi = x[i]; 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar yi = y[i]; 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x[i] = j.c() * xi + conj(j.s()) * yi; 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath y[i] = -j.s() * xi + conj(j.c()) * yi; 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /*** fixed-size vectorized path ***/ 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(VectorX::SizeAtCompileTime != Dynamic && 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath (VectorX::Flags & VectorY::Flags & PacketAccessBit) && 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath (VectorX::Flags & VectorY::Flags & AlignedBit)) 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Packet pc = pset1<Packet>(j.c()); 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Packet ps = pset1<Packet>(j.s()); 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath conj_helper<Packet,Packet,NumTraits<Scalar>::IsComplex,false> pcj; 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* EIGEN_RESTRICT px = x; 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* EIGEN_RESTRICT py = y; 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index i=0; i<size; i+=PacketSize) 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Packet xi = pload<Packet>(px); 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Packet yi = pload<Packet>(py); 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pstore(px, padd(pmul(pc,xi),pcj.pmul(ps,yi))); 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pstore(py, psub(pcj.pmul(pc,yi),pmul(ps,xi))); 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath px += PacketSize; 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath py += PacketSize; 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /*** non-vectorized path ***/ 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index i=0; i<size; ++i) 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar xi = *x; 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar yi = *y; 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *x = j.c() * xi + conj(j.s()) * yi; 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *y = -j.s() * xi + conj(j.c()) * yi; 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x += incrx; 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath y += incry; 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_JACOBI_H 421