Rev | Line | |
---|
[136] | 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>
|
---|
| 8 | void 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
|
---|
Note:
See
TracBrowser
for help on using the repository browser.