| namespace Eigen { |
| |
| namespace internal { |
| |
| template <typename Scalar> |
| void rwupdt( |
| Matrix< Scalar, Dynamic, Dynamic > &r, |
| const Matrix< Scalar, Dynamic, 1> &w, |
| Matrix< Scalar, Dynamic, 1> &b, |
| Scalar alpha) |
| { |
| typedef DenseIndex Index; |
| |
| const Index n = r.cols(); |
| eigen_assert(r.rows()>=n); |
| std::vector<JacobiRotation<Scalar> > givens(n); |
| |
| /* Local variables */ |
| Scalar temp, rowj; |
| |
| /* Function Body */ |
| for (Index j = 0; j < n; ++j) { |
| rowj = w[j]; |
| |
| /* apply the previous transformations to */ |
| /* r(i,j), i=0,1,...,j-1, and to w(j). */ |
| for (Index i = 0; i < j; ++i) { |
| temp = givens[i].c() * r(i,j) + givens[i].s() * rowj; |
| rowj = -givens[i].s() * r(i,j) + givens[i].c() * rowj; |
| r(i,j) = temp; |
| } |
| |
| /* determine a givens rotation which eliminates w(j). */ |
| givens[j].makeGivens(-r(j,j), rowj); |
| |
| if (rowj == 0.) |
| continue; // givens[j] is identity |
| |
| /* apply the current transformation to r(j,j), b(j), and alpha. */ |
| r(j,j) = givens[j].c() * r(j,j) + givens[j].s() * rowj; |
| temp = givens[j].c() * b[j] + givens[j].s() * alpha; |
| alpha = -givens[j].s() * b[j] + givens[j].c() * alpha; |
| b[j] = temp; |
| } |
| } |
| |
| } // end namespace internal |
| |
| } // end namespace Eigen |
