1 namespace Eigen {
2
3 namespace internal {
4
5 template <typename Scalar>
rwupdt(Matrix<Scalar,Dynamic,Dynamic> & r,const Matrix<Scalar,Dynamic,1> & w,Matrix<Scalar,Dynamic,1> & b,Scalar alpha)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