26 #ifndef EIGEN_JACOBI_H
27 #define EIGEN_JACOBI_H
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; }
78 template<
typename Derived>
82 void makeGivens(
const Scalar&
p,
const Scalar&
q, Scalar* z=0);
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);
96 template<
typename Scalar>
111 if(tau>RealScalar(0))
113 t = RealScalar(1) / (tau + w);
117 t = RealScalar(1) / (tau - w);
119 RealScalar sign_t = t > RealScalar(0) ? RealScalar(1) : RealScalar(-1);
136 template<
typename Scalar>
137 template<
typename Derived>
159 template<
typename Scalar>
167 template<
typename Scalar>
176 else if(p==Scalar(0))
222 template<
typename Scalar>
228 m_c = p<Scalar(0) ? Scalar(-1) : Scalar(1);
232 else if(p==Scalar(0))
235 m_s = q<Scalar(0) ? Scalar(1) : Scalar(-1);
272 template<
typename VectorX,
typename VectorY,
typename OtherScalar>
282 template<
typename Derived>
283 template<
typename OtherScalar>
286 RowXpr x(this->
row(p));
287 RowXpr
y(this->
row(q));
297 template<
typename Derived>
298 template<
typename OtherScalar>
301 ColXpr x(this->
col(p));
302 ColXpr
y(this->
col(q));
307 template<
typename VectorX,
typename VectorY,
typename OtherScalar>
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;
315 Index size = _x.size();
316 Index incrx = _x.innerStride();
317 Index incry = _y.innerStride();
324 if(VectorX::SizeAtCompileTime ==
Dynamic &&
326 ((incrx==1 && incry==1) || PacketSize == 1))
329 enum { Peeling = 2 };
331 Index alignedStart = internal::first_aligned(y, size);
332 Index alignedEnd = alignedStart + ((size-alignedStart)/PacketSize)*PacketSize;
334 const Packet pc = pset1<Packet>(j.
c());
335 const Packet ps = pset1<Packet>(j.
s());
338 for(Index i=0; i<alignedStart; ++i)
342 x[i] = j.
c() * xi +
conj(j.
s()) * yi;
343 y[i] = -j.
s() * xi +
conj(j.
c()) * yi;
349 if(internal::first_aligned(x, size)==alignedStart)
351 for(Index i=alignedStart; i<alignedEnd; i+=PacketSize)
353 Packet xi = pload<Packet>(px);
354 Packet yi = pload<Packet>(py);
363 Index peelingEnd = alignedStart + ((size-alignedStart)/(Peeling*PacketSize))*(Peeling*PacketSize);
364 for(Index i=alignedStart; i<peelingEnd; i+=Peeling*PacketSize)
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);
374 px += Peeling*PacketSize;
375 py += Peeling*PacketSize;
377 if(alignedEnd!=peelingEnd)
379 Packet xi = ploadu<Packet>(x+peelingEnd);
380 Packet yi = pload <Packet>(y+peelingEnd);
386 for(Index i=alignedEnd; i<size; ++i)
390 x[i] = j.
c() * xi +
conj(j.
s()) * yi;
391 y[i] = -j.
s() * xi +
conj(j.
c()) * yi;
396 else if(VectorX::SizeAtCompileTime !=
Dynamic &&
397 (VectorX::Flags & VectorY::Flags & PacketAccessBit) &&
398 (VectorX::Flags & VectorY::Flags &
AlignedBit))
400 const Packet pc = pset1<Packet>(j.
c());
401 const Packet ps = pset1<Packet>(j.
s());
405 for(Index i=0; i<size; i+=PacketSize)
407 Packet xi = pload<Packet>(px);
408 Packet yi = pload<Packet>(py);
419 for(Index i=0; i<size; ++i)
423 *x = j.
c() * xi +
conj(j.
s()) * yi;
424 *y = -j.
s() * xi +
conj(j.
c()) * yi;
435 #endif // EIGEN_JACOBI_H