1namespace Eigen {
2
3namespace internal {
4
5template <typename Scalar>
6void 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
100