You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
30 lines
1.1 KiB
30 lines
1.1 KiB
namespace Eigen {
|
|
|
|
namespace internal {
|
|
|
|
// TODO : move this to GivensQR once there's such a thing in Eigen
|
|
|
|
template <typename Scalar>
|
|
void r1mpyq(DenseIndex m, DenseIndex n, Scalar *a, const std::vector<JacobiRotation<Scalar> > &v_givens, const std::vector<JacobiRotation<Scalar> > &w_givens)
|
|
{
|
|
typedef DenseIndex Index;
|
|
|
|
/* apply the first set of givens rotations to a. */
|
|
for (Index j = n-2; j>=0; --j)
|
|
for (Index i = 0; i<m; ++i) {
|
|
Scalar temp = v_givens[j].c() * a[i+m*j] - v_givens[j].s() * a[i+m*(n-1)];
|
|
a[i+m*(n-1)] = v_givens[j].s() * a[i+m*j] + v_givens[j].c() * a[i+m*(n-1)];
|
|
a[i+m*j] = temp;
|
|
}
|
|
/* apply the second set of givens rotations to a. */
|
|
for (Index j = 0; j<n-1; ++j)
|
|
for (Index i = 0; i<m; ++i) {
|
|
Scalar temp = w_givens[j].c() * a[i+m*j] + w_givens[j].s() * a[i+m*(n-1)];
|
|
a[i+m*(n-1)] = -w_givens[j].s() * a[i+m*j] + w_givens[j].c() * a[i+m*(n-1)];
|
|
a[i+m*j] = temp;
|
|
}
|
|
}
|
|
|
|
} // end namespace internal
|
|
|
|
} // end namespace Eigen
|
|
|