Jacobi.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
5 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
6 //
7 // Eigen is free software; you can redistribute it and/or
8 // modify it under the terms of the GNU Lesser General Public
9 // License as published by the Free Software Foundation; either
10 // version 3 of the License, or (at your option) any later version.
11 //
12 // Alternatively, you can redistribute it and/or
13 // modify it under the terms of the GNU General Public License as
14 // published by the Free Software Foundation; either version 2 of
15 // the License, or (at your option) any later version.
16 //
17 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
18 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
20 // GNU General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public
23 // License and a copy of the GNU General Public License along with
24 // Eigen. If not, see <http://www.gnu.org/licenses/>.
25 
26 #ifndef EIGEN_JACOBI_H
27 #define EIGEN_JACOBI_H
28 
29 namespace Eigen {
30 
49 template<typename Scalar> class JacobiRotation
50 {
51  public:
53 
56 
58  JacobiRotation(const Scalar& c, const Scalar& s) : m_c(c), m_s(s) {}
59 
60  Scalar& c() { return m_c; }
61  Scalar c() const { return m_c; }
62  Scalar& s() { return m_s; }
63  Scalar s() const { return m_s; }
64 
67  {
68  return JacobiRotation(m_c * other.m_c - internal::conj(m_s) * other.m_s,
70  }
71 
74 
77 
78  template<typename Derived>
79  bool makeJacobi(const MatrixBase<Derived>&, typename Derived::Index p, typename Derived::Index q);
80  bool makeJacobi(RealScalar x, Scalar y, RealScalar z);
81 
82  void makeGivens(const Scalar& p, const Scalar& q, Scalar* z=0);
83 
84  protected:
85  void makeGivens(const Scalar& p, const Scalar& q, Scalar* z, internal::true_type);
86  void makeGivens(const Scalar& p, const Scalar& q, Scalar* z, internal::false_type);
87 
88  Scalar m_c, m_s;
89 };
90 
96 template<typename Scalar>
98 {
99  typedef typename NumTraits<Scalar>::Real RealScalar;
100  if(y == Scalar(0))
101  {
102  m_c = Scalar(1);
103  m_s = Scalar(0);
104  return false;
105  }
106  else
107  {
108  RealScalar tau = (x-z)/(RealScalar(2)*internal::abs(y));
109  RealScalar w = internal::sqrt(internal::abs2(tau) + RealScalar(1));
110  RealScalar t;
111  if(tau>RealScalar(0))
112  {
113  t = RealScalar(1) / (tau + w);
114  }
115  else
116  {
117  t = RealScalar(1) / (tau - w);
118  }
119  RealScalar sign_t = t > RealScalar(0) ? RealScalar(1) : RealScalar(-1);
120  RealScalar n = RealScalar(1) / internal::sqrt(internal::abs2(t)+RealScalar(1));
121  m_s = - sign_t * (internal::conj(y) / internal::abs(y)) * internal::abs(t) * n;
122  m_c = n;
123  return true;
124  }
125 }
126 
136 template<typename Scalar>
137 template<typename Derived>
138 inline bool JacobiRotation<Scalar>::makeJacobi(const MatrixBase<Derived>& m, typename Derived::Index p, typename Derived::Index q)
139 {
140  return makeJacobi(internal::real(m.coeff(p,p)), m.coeff(p,q), internal::real(m.coeff(q,q)));
141 }
142 
159 template<typename Scalar>
160 void JacobiRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* z)
161 {
162  makeGivens(p, q, z, typename internal::conditional<NumTraits<Scalar>::IsComplex, internal::true_type, internal::false_type>::type());
163 }
164 
165 
166 // specialization for complexes
167 template<typename Scalar>
168 void JacobiRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* r, internal::true_type)
169 {
170  if(q==Scalar(0))
171  {
172  m_c = internal::real(p)<0 ? Scalar(-1) : Scalar(1);
173  m_s = 0;
174  if(r) *r = m_c * p;
175  }
176  else if(p==Scalar(0))
177  {
178  m_c = 0;
179  m_s = -q/internal::abs(q);
180  if(r) *r = internal::abs(q);
181  }
182  else
183  {
184  RealScalar p1 = internal::norm1(p);
185  RealScalar q1 = internal::norm1(q);
186  if(p1>=q1)
187  {
188  Scalar ps = p / p1;
189  RealScalar p2 = internal::abs2(ps);
190  Scalar qs = q / p1;
191  RealScalar q2 = internal::abs2(qs);
192 
193  RealScalar u = internal::sqrt(RealScalar(1) + q2/p2);
194  if(internal::real(p)<RealScalar(0))
195  u = -u;
196 
197  m_c = Scalar(1)/u;
198  m_s = -qs*internal::conj(ps)*(m_c/p2);
199  if(r) *r = p * u;
200  }
201  else
202  {
203  Scalar ps = p / q1;
204  RealScalar p2 = internal::abs2(ps);
205  Scalar qs = q / q1;
206  RealScalar q2 = internal::abs2(qs);
207 
208  RealScalar u = q1 * internal::sqrt(p2 + q2);
209  if(internal::real(p)<RealScalar(0))
210  u = -u;
211 
212  p1 = internal::abs(p);
213  ps = p/p1;
214  m_c = p1/u;
215  m_s = -internal::conj(ps) * (q/u);
216  if(r) *r = ps * u;
217  }
218  }
219 }
220 
221 // specialization for reals
222 template<typename Scalar>
223 void JacobiRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* r, internal::false_type)
224 {
225 
226  if(q==Scalar(0))
227  {
228  m_c = p<Scalar(0) ? Scalar(-1) : Scalar(1);
229  m_s = Scalar(0);
230  if(r) *r = internal::abs(p);
231  }
232  else if(p==Scalar(0))
233  {
234  m_c = Scalar(0);
235  m_s = q<Scalar(0) ? Scalar(1) : Scalar(-1);
236  if(r) *r = internal::abs(q);
237  }
238  else if(internal::abs(p) > internal::abs(q))
239  {
240  Scalar t = q/p;
241  Scalar u = internal::sqrt(Scalar(1) + internal::abs2(t));
242  if(p<Scalar(0))
243  u = -u;
244  m_c = Scalar(1)/u;
245  m_s = -t * m_c;
246  if(r) *r = p * u;
247  }
248  else
249  {
250  Scalar t = p/q;
251  Scalar u = internal::sqrt(Scalar(1) + internal::abs2(t));
252  if(q<Scalar(0))
253  u = -u;
254  m_s = -Scalar(1)/u;
255  m_c = -t * m_s;
256  if(r) *r = q * u;
257  }
258 
259 }
260 
261 /****************************************************************************************
262 * Implementation of MatrixBase methods
263 ****************************************************************************************/
264 
271 namespace internal {
272 template<typename VectorX, typename VectorY, typename OtherScalar>
273 void apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, const JacobiRotation<OtherScalar>& j);
274 }
275 
282 template<typename Derived>
283 template<typename OtherScalar>
285 {
286  RowXpr x(this->row(p));
287  RowXpr y(this->row(q));
289 }
290 
297 template<typename Derived>
298 template<typename OtherScalar>
300 {
301  ColXpr x(this->col(p));
302  ColXpr y(this->col(q));
304 }
305 
306 namespace internal {
307 template<typename VectorX, typename VectorY, typename OtherScalar>
308 void /*EIGEN_DONT_INLINE*/ apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, const JacobiRotation<OtherScalar>& j)
309 {
310  typedef typename VectorX::Index Index;
311  typedef typename VectorX::Scalar Scalar;
312  enum { PacketSize = packet_traits<Scalar>::size };
313  typedef typename packet_traits<Scalar>::type Packet;
314  eigen_assert(_x.size() == _y.size());
315  Index size = _x.size();
316  Index incrx = _x.innerStride();
317  Index incry = _y.innerStride();
318 
319  Scalar* EIGEN_RESTRICT x = &_x.coeffRef(0);
320  Scalar* EIGEN_RESTRICT y = &_y.coeffRef(0);
321 
322  /*** dynamic-size vectorized paths ***/
323 
324  if(VectorX::SizeAtCompileTime == Dynamic &&
325  (VectorX::Flags & VectorY::Flags & PacketAccessBit) &&
326  ((incrx==1 && incry==1) || PacketSize == 1))
327  {
328  // both vectors are sequentially stored in memory => vectorization
329  enum { Peeling = 2 };
330 
331  Index alignedStart = internal::first_aligned(y, size);
332  Index alignedEnd = alignedStart + ((size-alignedStart)/PacketSize)*PacketSize;
333 
334  const Packet pc = pset1<Packet>(j.c());
335  const Packet ps = pset1<Packet>(j.s());
337 
338  for(Index i=0; i<alignedStart; ++i)
339  {
340  Scalar xi = x[i];
341  Scalar yi = y[i];
342  x[i] = j.c() * xi + conj(j.s()) * yi;
343  y[i] = -j.s() * xi + conj(j.c()) * yi;
344  }
345 
346  Scalar* EIGEN_RESTRICT px = x + alignedStart;
347  Scalar* EIGEN_RESTRICT py = y + alignedStart;
348 
349  if(internal::first_aligned(x, size)==alignedStart)
350  {
351  for(Index i=alignedStart; i<alignedEnd; i+=PacketSize)
352  {
353  Packet xi = pload<Packet>(px);
354  Packet yi = pload<Packet>(py);
355  pstore(px, padd(pmul(pc,xi),pcj.pmul(ps,yi)));
356  pstore(py, psub(pcj.pmul(pc,yi),pmul(ps,xi)));
357  px += PacketSize;
358  py += PacketSize;
359  }
360  }
361  else
362  {
363  Index peelingEnd = alignedStart + ((size-alignedStart)/(Peeling*PacketSize))*(Peeling*PacketSize);
364  for(Index i=alignedStart; i<peelingEnd; i+=Peeling*PacketSize)
365  {
366  Packet xi = ploadu<Packet>(px);
367  Packet xi1 = ploadu<Packet>(px+PacketSize);
368  Packet yi = pload <Packet>(py);
369  Packet yi1 = pload <Packet>(py+PacketSize);
370  pstoreu(px, padd(pmul(pc,xi),pcj.pmul(ps,yi)));
371  pstoreu(px+PacketSize, padd(pmul(pc,xi1),pcj.pmul(ps,yi1)));
372  pstore (py, psub(pcj.pmul(pc,yi),pmul(ps,xi)));
373  pstore (py+PacketSize, psub(pcj.pmul(pc,yi1),pmul(ps,xi1)));
374  px += Peeling*PacketSize;
375  py += Peeling*PacketSize;
376  }
377  if(alignedEnd!=peelingEnd)
378  {
379  Packet xi = ploadu<Packet>(x+peelingEnd);
380  Packet yi = pload <Packet>(y+peelingEnd);
381  pstoreu(x+peelingEnd, padd(pmul(pc,xi),pcj.pmul(ps,yi)));
382  pstore (y+peelingEnd, psub(pcj.pmul(pc,yi),pmul(ps,xi)));
383  }
384  }
385 
386  for(Index i=alignedEnd; i<size; ++i)
387  {
388  Scalar xi = x[i];
389  Scalar yi = y[i];
390  x[i] = j.c() * xi + conj(j.s()) * yi;
391  y[i] = -j.s() * xi + conj(j.c()) * yi;
392  }
393  }
394 
395  /*** fixed-size vectorized path ***/
396  else if(VectorX::SizeAtCompileTime != Dynamic &&
397  (VectorX::Flags & VectorY::Flags & PacketAccessBit) &&
398  (VectorX::Flags & VectorY::Flags & AlignedBit))
399  {
400  const Packet pc = pset1<Packet>(j.c());
401  const Packet ps = pset1<Packet>(j.s());
403  Scalar* EIGEN_RESTRICT px = x;
404  Scalar* EIGEN_RESTRICT py = y;
405  for(Index i=0; i<size; i+=PacketSize)
406  {
407  Packet xi = pload<Packet>(px);
408  Packet yi = pload<Packet>(py);
409  pstore(px, padd(pmul(pc,xi),pcj.pmul(ps,yi)));
410  pstore(py, psub(pcj.pmul(pc,yi),pmul(ps,xi)));
411  px += PacketSize;
412  py += PacketSize;
413  }
414  }
415 
416  /*** non-vectorized path ***/
417  else
418  {
419  for(Index i=0; i<size; ++i)
420  {
421  Scalar xi = *x;
422  Scalar yi = *y;
423  *x = j.c() * xi + conj(j.s()) * yi;
424  *y = -j.s() * xi + conj(j.c()) * yi;
425  x += incrx;
426  y += incry;
427  }
428  }
429 }
430 
431 } // end namespace internal
432 
433 } // end namespace Eigen
434 
435 #endif // EIGEN_JACOBI_H