1 namespace Eigen { 2 3 namespace internal { 4 5 template <typename Scalar> 6 void rwupdt( 7 Matrix< Scalar, Dynamic, Dynamic > &r, 8 const Matrix< Scalar, Dynamic, 1> &w, 9 Matrix< Scalar, Dynamic, 1> &b, 10 Scalar alpha) 11 { 12 typedef DenseIndex Index; 13 14 const Index n = r.cols(); 15 eigen_assert(r.rows()>=n); 16 std::vector<JacobiRotation<Scalar> > givens(n); 17 18 /* Local variables */ 19 Scalar temp, rowj; 20 21 /* Function Body */ 22 for (Index j = 0; j < n; ++j) { 23 rowj = w[j]; 24 25 /* apply the previous transformations to */ 26 /* r(i,j), i=0,1,...,j-1, and to w(j). */ 27 for (Index i = 0; i < j; ++i) { 28 temp = givens[i].c() * r(i,j) + givens[i].s() * rowj; 29 rowj = -givens[i].s() * r(i,j) + givens[i].c() * rowj; 30 r(i,j) = temp; 31 } 32 33 /* determine a givens rotation which eliminates w(j). */ 34 givens[j].makeGivens(-r(j,j), rowj); 35 36 if (rowj == 0.) 37 continue; // givens[j] is identity 38 39 /* apply the current transformation to r(j,j), b(j), and alpha. */ 40 r(j,j) = givens[j].c() * r(j,j) + givens[j].s() * rowj; 41 temp = givens[j].c() * b[j] + givens[j].s() * alpha; 42 alpha = -givens[j].s() * b[j] + givens[j].c() * alpha; 43 b[j] = temp; 44 } 45 } 46 47 } // end namespace internal 48 49 } // end namespace Eigen 50
