src/NonLinearOptimization/r1updt.h

0001 namespace Eigen {
0002
0003 namespace internal {
0004
0005 template <typename Scalar>
0006 void r1updt(
0007         Matrix< Scalar, Dynamic, Dynamic > &s,
0008         const Matrix< Scalar, Dynamic, 1> &u,
0009         std::vector<JacobiRotation<Scalar> > &v_givens,
0010         std::vector<JacobiRotation<Scalar> > &w_givens,
0011         Matrix< Scalar, Dynamic, 1> &v,
0012         Matrix< Scalar, Dynamic, 1> &w,
0013         bool *sing)
0014 {
0015     typedef DenseIndex Index;
0016     const JacobiRotation<Scalar> IdentityRotation = JacobiRotation<Scalar>(1,0);
0017
0018     /* Local variables */
0019     const Index m = s.rows();
0020     const Index n = s.cols();
0021     Index i, j=1;
0022     Scalar temp;
0023     JacobiRotation<Scalar> givens;
0024
0025     // r1updt had a broader usecase, but we don't use it here. And, more
0026     // importantly, we can not test it.
0027     eigen_assert(m==n);
0028     eigen_assert(u.size()==m);
0029     eigen_assert(v.size()==n);
0030     eigen_assert(w.size()==n);
0031
0032     /* move the nontrivial part of the last column of s into w. */
0033     w[n-1] = s(n-1,n-1);
0034
0035     /* rotate the vector v into a multiple of the n-th unit vector */
0036     /* in such a way that a spike is introduced into w. */
0037     for (j=n-2; j>=0; --j) {
0038         w[j] = 0.;
0039         if (v[j] != 0.) {
0040             /* determine a givens rotation which eliminates the */
0041             /* j-th element of v. */
0042             givens.makeGivens(-v[n-1], v[j]);
0043
0044             /* apply the transformation to v and store the information */
0045             /* necessary to recover the givens rotation. */
0046             v[n-1] = givens.s() * v[j] + givens.c() * v[n-1];
0047             v_givens[j] = givens;
0048
0049             /* apply the transformation to s and extend the spike in w. */
0050             for (i = j; i < m; ++i) {
0051                 temp = givens.c() * s(j,i) - givens.s() * w[i];
0052                 w[i] = givens.s() * s(j,i) + givens.c() * w[i];
0053                 s(j,i) = temp;
0054             }
0055         } else
0056             v_givens[j] = IdentityRotation;
0057     }
0058
0059     /* add the spike from the rank 1 update to w. */
0060     w += v[n-1] * u;
0061
0062     /* eliminate the spike. */
0063     *sing = false;
0064     for (j = 0; j < n-1; ++j) {
0065         if (w[j] != 0.) {
0066             /* determine a givens rotation which eliminates the */
0067             /* j-th element of the spike. */
0068             givens.makeGivens(-s(j,j), w[j]);
0069
0070             /* apply the transformation to s and reduce the spike in w. */
0071             for (i = j; i < m; ++i) {
0072                 temp = givens.c() * s(j,i) + givens.s() * w[i];
0073                 w[i] = -givens.s() * s(j,i) + givens.c() * w[i];
0074                 s(j,i) = temp;
0075             }
0076
0077             /* store the information necessary to recover the */
0078             /* givens rotation. */
0079             w_givens[j] = givens;
0080         } else
0081             v_givens[j] = IdentityRotation;
0082
0083         /* test for zero diagonal elements in the output s. */
0084         if (s(j,j) == 0.) {
0085             *sing = true;
0086         }
0087     }
0088     /* move w back into the last column of the output s. */
0089     s(n-1,n-1) = w[n-1];
0090
0091     if (s(j,j) == 0.) {
0092         *sing = true;
0093     }
0094     return;
0095 }
0096
0097 } // end namespace internal
0098
0099 } // end namespace Eigen