1 namespace Eigen {
2
3 namespace internal {
4
5 template <typename Scalar>
r1updt(Matrix<Scalar,Dynamic,Dynamic> & s,const Matrix<Scalar,Dynamic,1> & u,std::vector<JacobiRotation<Scalar>> & v_givens,std::vector<JacobiRotation<Scalar>> & w_givens,Matrix<Scalar,Dynamic,1> & v,Matrix<Scalar,Dynamic,1> & w,bool * sing)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
100