| 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 |