[136] | 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
|
---|