1namespace Eigen {
2
3namespace internal {
4
5// TODO : move this to GivensQR once there's such a thing in Eigen
6
7template <typename Scalar>
8void r1mpyq(DenseIndex m, DenseIndex n, Scalar *a, const std::vector<JacobiRotation<Scalar> > &v_givens, const std::vector<JacobiRotation<Scalar> > &w_givens)
9{
10    typedef DenseIndex Index;
11
12    /*     apply the first set of givens rotations to a. */
13    for (Index j = n-2; j>=0; --j)
14        for (Index i = 0; i<m; ++i) {
15            Scalar temp = v_givens[j].c() * a[i+m*j] - v_givens[j].s() * a[i+m*(n-1)];
16            a[i+m*(n-1)] = v_givens[j].s() * a[i+m*j] + v_givens[j].c() * a[i+m*(n-1)];
17            a[i+m*j] = temp;
18        }
19    /*     apply the second set of givens rotations to a. */
20    for (Index j = 0; j<n-1; ++j)
21        for (Index i = 0; i<m; ++i) {
22            Scalar temp = w_givens[j].c() * a[i+m*j] + w_givens[j].s() * a[i+m*(n-1)];
23            a[i+m*(n-1)] = -w_givens[j].s() * a[i+m*j] + w_givens[j].c() * a[i+m*(n-1)];
24            a[i+m*j] = temp;
25        }
26}
27
28} // end namespace internal
29
30} // end namespace Eigen
31