| 1 | namespace Eigen {
|
|---|
| 2 |
|
|---|
| 3 | namespace internal {
|
|---|
| 4 |
|
|---|
| 5 | template <typename Scalar>
|
|---|
| 6 | void r1updt(
|
|---|
| 7 | Matrix< Scalar, Dynamic, Dynamic > &s,
|
|---|
| 8 | const Matrix< Scalar, Dynamic, 1> &u,
|
|---|
| 9 | std::vector<JacobiRotation<Scalar> > &v_givens,
|
|---|
| 10 | std::vector<JacobiRotation<Scalar> > &w_givens,
|
|---|
| 11 | Matrix< Scalar, Dynamic, 1> &v,
|
|---|
| 12 | Matrix< Scalar, Dynamic, 1> &w,
|
|---|
| 13 | bool *sing)
|
|---|
| 14 | {
|
|---|
| 15 | typedef DenseIndex Index;
|
|---|
| 16 | const JacobiRotation<Scalar> IdentityRotation = JacobiRotation<Scalar>(1,0);
|
|---|
| 17 |
|
|---|
| 18 | /* Local variables */
|
|---|
| 19 | const Index m = s.rows();
|
|---|
| 20 | const Index n = s.cols();
|
|---|
| 21 | Index i, j=1;
|
|---|
| 22 | Scalar temp;
|
|---|
| 23 | JacobiRotation<Scalar> givens;
|
|---|
| 24 |
|
|---|
| 25 | // r1updt had a broader usecase, but we dont use it here. And, more
|
|---|
| 26 | // importantly, we can not test it.
|
|---|
| 27 | eigen_assert(m==n);
|
|---|
| 28 | eigen_assert(u.size()==m);
|
|---|
| 29 | eigen_assert(v.size()==n);
|
|---|
| 30 | eigen_assert(w.size()==n);
|
|---|
| 31 |
|
|---|
| 32 | /* move the nontrivial part of the last column of s into w. */
|
|---|
| 33 | w[n-1] = s(n-1,n-1);
|
|---|
| 34 |
|
|---|
| 35 | /* rotate the vector v into a multiple of the n-th unit vector */
|
|---|
| 36 | /* in such a way that a spike is introduced into w. */
|
|---|
| 37 | for (j=n-2; j>=0; --j) {
|
|---|
| 38 | w[j] = 0.;
|
|---|
| 39 | if (v[j] != 0.) {
|
|---|
| 40 | /* determine a givens rotation which eliminates the */
|
|---|
| 41 | /* j-th element of v. */
|
|---|
| 42 | givens.makeGivens(-v[n-1], v[j]);
|
|---|
| 43 |
|
|---|
| 44 | /* apply the transformation to v and store the information */
|
|---|
| 45 | /* necessary to recover the givens rotation. */
|
|---|
| 46 | v[n-1] = givens.s() * v[j] + givens.c() * v[n-1];
|
|---|
| 47 | v_givens[j] = givens;
|
|---|
| 48 |
|
|---|
| 49 | /* apply the transformation to s and extend the spike in w. */
|
|---|
| 50 | for (i = j; i < m; ++i) {
|
|---|
| 51 | temp = givens.c() * s(j,i) - givens.s() * w[i];
|
|---|
| 52 | w[i] = givens.s() * s(j,i) + givens.c() * w[i];
|
|---|
| 53 | s(j,i) = temp;
|
|---|
| 54 | }
|
|---|
| 55 | } else
|
|---|
| 56 | v_givens[j] = IdentityRotation;
|
|---|
| 57 | }
|
|---|
| 58 |
|
|---|
| 59 | /* add the spike from the rank 1 update to w. */
|
|---|
| 60 | w += v[n-1] * u;
|
|---|
| 61 |
|
|---|
| 62 | /* eliminate the spike. */
|
|---|
| 63 | *sing = false;
|
|---|
| 64 | for (j = 0; j < n-1; ++j) {
|
|---|
| 65 | if (w[j] != 0.) {
|
|---|
| 66 | /* determine a givens rotation which eliminates the */
|
|---|
| 67 | /* j-th element of the spike. */
|
|---|
| 68 | givens.makeGivens(-s(j,j), w[j]);
|
|---|
| 69 |
|
|---|
| 70 | /* apply the transformation to s and reduce the spike in w. */
|
|---|
| 71 | for (i = j; i < m; ++i) {
|
|---|
| 72 | temp = givens.c() * s(j,i) + givens.s() * w[i];
|
|---|
| 73 | w[i] = -givens.s() * s(j,i) + givens.c() * w[i];
|
|---|
| 74 | s(j,i) = temp;
|
|---|
| 75 | }
|
|---|
| 76 |
|
|---|
| 77 | /* store the information necessary to recover the */
|
|---|
| 78 | /* givens rotation. */
|
|---|
| 79 | w_givens[j] = givens;
|
|---|
| 80 | } else
|
|---|
| 81 | v_givens[j] = IdentityRotation;
|
|---|
| 82 |
|
|---|
| 83 | /* test for zero diagonal elements in the output s. */
|
|---|
| 84 | if (s(j,j) == 0.) {
|
|---|
| 85 | *sing = true;
|
|---|
| 86 | }
|
|---|
| 87 | }
|
|---|
| 88 | /* move w back into the last column of the output s. */
|
|---|
| 89 | s(n-1,n-1) = w[n-1];
|
|---|
| 90 |
|
|---|
| 91 | if (s(j,j) == 0.) {
|
|---|
| 92 | *sing = true;
|
|---|
| 93 | }
|
|---|
| 94 | return;
|
|---|
| 95 | }
|
|---|
| 96 |
|
|---|
| 97 | } // end namespace internal
|
|---|
| 98 |
|
|---|
| 99 | } // end namespace Eigen
|
|---|
