11 #ifndef EIGEN_JACOBI_H
12 #define EIGEN_JACOBI_H
45 Scalar& c() {
return m_c; }
46 Scalar c()
const {
return m_c; }
47 Scalar& s() {
return m_s; }
48 Scalar s()
const {
return m_s; }
53 return JacobiRotation(m_c * other.m_c - internal::conj(m_s) * other.m_s,
54 internal::conj(m_c * internal::conj(other.m_s) + internal::conj(m_s) * internal::conj(other.m_c)));
63 template<
typename Derived>
65 bool makeJacobi(RealScalar x, Scalar y, RealScalar z);
67 void makeGivens(
const Scalar& p,
const Scalar& q, Scalar* z=0);
70 void makeGivens(
const Scalar& p,
const Scalar& q, Scalar* z, internal::true_type);
71 void makeGivens(
const Scalar& p,
const Scalar& q, Scalar* z, internal::false_type);
81 template<
typename Scalar>
93 RealScalar tau = (x-z)/(RealScalar(2)*internal::abs(y));
94 RealScalar w = internal::sqrt(internal::abs2(tau) + RealScalar(1));
98 t = RealScalar(1) / (tau + w);
102 t = RealScalar(1) / (tau - w);
104 RealScalar sign_t = t > RealScalar(0) ? RealScalar(1) : RealScalar(-1);
105 RealScalar n = RealScalar(1) / internal::sqrt(internal::abs2(t)+RealScalar(1));
106 m_s = - sign_t * (internal::conj(y) / internal::abs(y)) * internal::abs(t) * n;
121 template<
typename Scalar>
122 template<
typename Derived>
125 return makeJacobi(internal::real(m.coeff(p,p)), m.coeff(p,q), internal::real(m.coeff(q,q)));
144 template<
typename Scalar>
152 template<
typename Scalar>
157 m_c = internal::real(p)<0 ? Scalar(-1) : Scalar(1);
161 else if(p==Scalar(0))
164 m_s = -q/internal::abs(q);
165 if(r) *r = internal::abs(q);
169 RealScalar p1 = internal::norm1(p);
170 RealScalar q1 = internal::norm1(q);
174 RealScalar p2 = internal::abs2(ps);
176 RealScalar q2 = internal::abs2(qs);
178 RealScalar u = internal::sqrt(RealScalar(1) + q2/p2);
179 if(internal::real(p)<RealScalar(0))
183 m_s = -qs*internal::conj(ps)*(m_c/p2);
189 RealScalar p2 = internal::abs2(ps);
191 RealScalar q2 = internal::abs2(qs);
193 RealScalar u = q1 * internal::sqrt(p2 + q2);
194 if(internal::real(p)<RealScalar(0))
197 p1 = internal::abs(p);
200 m_s = -internal::conj(ps) * (q/u);
207 template<
typename Scalar>
213 m_c = p<Scalar(0) ? Scalar(-1) : Scalar(1);
215 if(r) *r = internal::abs(p);
217 else if(p==Scalar(0))
220 m_s = q<Scalar(0) ? Scalar(1) : Scalar(-1);
221 if(r) *r = internal::abs(q);
223 else if(internal::abs(p) > internal::abs(q))
226 Scalar u = internal::sqrt(Scalar(1) + internal::abs2(t));
236 Scalar u = internal::sqrt(Scalar(1) + internal::abs2(t));
257 template<
typename VectorX,
typename VectorY,
typename OtherScalar>
258 void apply_rotation_in_the_plane(VectorX& _x, VectorY& _y,
const JacobiRotation<OtherScalar>& j);
267 template<
typename Derived>
268 template<
typename OtherScalar>
271 RowXpr x(this->row(p));
272 RowXpr y(this->row(q));
273 internal::apply_rotation_in_the_plane(x, y, j);
282 template<
typename Derived>
283 template<
typename OtherScalar>
286 ColXpr x(this->col(p));
287 ColXpr y(this->col(q));
288 internal::apply_rotation_in_the_plane(x, y, j.
transpose());
292 template<
typename VectorX,
typename VectorY,
typename OtherScalar>
295 typedef typename VectorX::Index Index;
296 typedef typename VectorX::Scalar Scalar;
297 enum { PacketSize = packet_traits<Scalar>::size };
298 typedef typename packet_traits<Scalar>::type Packet;
299 eigen_assert(_x.size() == _y.size());
300 Index size = _x.size();
301 Index incrx = _x.innerStride();
302 Index incry = _y.innerStride();
304 Scalar* EIGEN_RESTRICT x = &_x.coeffRef(0);
305 Scalar* EIGEN_RESTRICT y = &_y.coeffRef(0);
309 if(VectorX::SizeAtCompileTime ==
Dynamic &&
311 ((incrx==1 && incry==1) || PacketSize == 1))
314 enum { Peeling = 2 };
316 Index alignedStart = internal::first_aligned(y, size);
317 Index alignedEnd = alignedStart + ((size-alignedStart)/PacketSize)*PacketSize;
319 const Packet pc = pset1<Packet>(j.c());
320 const Packet ps = pset1<Packet>(j.s());
321 conj_helper<Packet,Packet,NumTraits<Scalar>::IsComplex,
false> pcj;
323 for(Index i=0; i<alignedStart; ++i)
327 x[i] = j.c() * xi + conj(j.s()) * yi;
328 y[i] = -j.s() * xi + conj(j.c()) * yi;
331 Scalar* EIGEN_RESTRICT px = x + alignedStart;
332 Scalar* EIGEN_RESTRICT py = y + alignedStart;
334 if(internal::first_aligned(x, size)==alignedStart)
336 for(Index i=alignedStart; i<alignedEnd; i+=PacketSize)
338 Packet xi = pload<Packet>(px);
339 Packet yi = pload<Packet>(py);
340 pstore(px, padd(pmul(pc,xi),pcj.pmul(ps,yi)));
341 pstore(py, psub(pcj.pmul(pc,yi),pmul(ps,xi)));
348 Index peelingEnd = alignedStart + ((size-alignedStart)/(Peeling*PacketSize))*(Peeling*PacketSize);
349 for(Index i=alignedStart; i<peelingEnd; i+=Peeling*PacketSize)
351 Packet xi = ploadu<Packet>(px);
352 Packet xi1 = ploadu<Packet>(px+PacketSize);
353 Packet yi = pload <Packet>(py);
354 Packet yi1 = pload <Packet>(py+PacketSize);
355 pstoreu(px, padd(pmul(pc,xi),pcj.pmul(ps,yi)));
356 pstoreu(px+PacketSize, padd(pmul(pc,xi1),pcj.pmul(ps,yi1)));
357 pstore (py, psub(pcj.pmul(pc,yi),pmul(ps,xi)));
358 pstore (py+PacketSize, psub(pcj.pmul(pc,yi1),pmul(ps,xi1)));
359 px += Peeling*PacketSize;
360 py += Peeling*PacketSize;
362 if(alignedEnd!=peelingEnd)
364 Packet xi = ploadu<Packet>(x+peelingEnd);
365 Packet yi = pload <Packet>(y+peelingEnd);
366 pstoreu(x+peelingEnd, padd(pmul(pc,xi),pcj.pmul(ps,yi)));
367 pstore (y+peelingEnd, psub(pcj.pmul(pc,yi),pmul(ps,xi)));
371 for(Index i=alignedEnd; i<size; ++i)
375 x[i] = j.c() * xi + conj(j.s()) * yi;
376 y[i] = -j.s() * xi + conj(j.c()) * yi;
381 else if(VectorX::SizeAtCompileTime !=
Dynamic &&
382 (VectorX::Flags & VectorY::Flags & PacketAccessBit) &&
383 (VectorX::Flags & VectorY::Flags &
AlignedBit))
385 const Packet pc = pset1<Packet>(j.c());
386 const Packet ps = pset1<Packet>(j.s());
387 conj_helper<Packet,Packet,NumTraits<Scalar>::IsComplex,
false> pcj;
388 Scalar* EIGEN_RESTRICT px = x;
389 Scalar* EIGEN_RESTRICT py = y;
390 for(Index i=0; i<size; i+=PacketSize)
392 Packet xi = pload<Packet>(px);
393 Packet yi = pload<Packet>(py);
394 pstore(px, padd(pmul(pc,xi),pcj.pmul(ps,yi)));
395 pstore(py, psub(pcj.pmul(pc,yi),pmul(ps,xi)));
404 for(Index i=0; i<size; ++i)
408 *x = j.c() * xi + conj(j.s()) * yi;
409 *y = -j.s() * xi + conj(j.c()) * yi;
420 #endif // EIGEN_JACOBI_H