1namespace Eigen {
2
3namespace internal {
4
5template <typename Scalar>
6void 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