1 namespace Eigen { 2 3 namespace internal { 4 5 // TODO : move this to GivensQR once there's such a thing in Eigen 6 7 template <typename Scalar> r1mpyq(DenseIndex m,DenseIndex n,Scalar * a,const std::vector<JacobiRotation<Scalar>> & v_givens,const std::vector<JacobiRotation<Scalar>> & w_givens)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