• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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