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 100
