r1updt.h
Go to the documentation of this file.
00001 namespace Eigen { 
00002 
00003 namespace internal {
00004 
00005 template <typename Scalar>
00006 void r1updt(
00007         Matrix< Scalar, Dynamic, Dynamic > &s,
00008         const Matrix< Scalar, Dynamic, 1> &u,
00009         std::vector<JacobiRotation<Scalar> > &v_givens,
00010         std::vector<JacobiRotation<Scalar> > &w_givens,
00011         Matrix< Scalar, Dynamic, 1> &v,
00012         Matrix< Scalar, Dynamic, 1> &w,
00013         bool *sing)
00014 {
00015     typedef DenseIndex Index;
00016     const JacobiRotation<Scalar> IdentityRotation = JacobiRotation<Scalar>(1,0);
00017 
00018     /* Local variables */
00019     const Index m = s.rows();
00020     const Index n = s.cols();
00021     Index i, j=1;
00022     Scalar temp;
00023     JacobiRotation<Scalar> givens;
00024 
00025     // r1updt had a broader usecase, but we dont use it here. And, more
00026     // importantly, we can not test it.
00027     eigen_assert(m==n);
00028     eigen_assert(u.size()==m);
00029     eigen_assert(v.size()==n);
00030     eigen_assert(w.size()==n);
00031 
00032     /* move the nontrivial part of the last column of s into w. */
00033     w[n-1] = s(n-1,n-1);
00034 
00035     /* rotate the vector v into a multiple of the n-th unit vector */
00036     /* in such a way that a spike is introduced into w. */
00037     for (j=n-2; j>=0; --j) {
00038         w[j] = 0.;
00039         if (v[j] != 0.) {
00040             /* determine a givens rotation which eliminates the */
00041             /* j-th element of v. */
00042             givens.makeGivens(-v[n-1], v[j]);
00043 
00044             /* apply the transformation to v and store the information */
00045             /* necessary to recover the givens rotation. */
00046             v[n-1] = givens.s() * v[j] + givens.c() * v[n-1];
00047             v_givens[j] = givens;
00048 
00049             /* apply the transformation to s and extend the spike in w. */
00050             for (i = j; i < m; ++i) {
00051                 temp = givens.c() * s(j,i) - givens.s() * w[i];
00052                 w[i] = givens.s() * s(j,i) + givens.c() * w[i];
00053                 s(j,i) = temp;
00054             }
00055         } else
00056             v_givens[j] = IdentityRotation;
00057     }
00058 
00059     /* add the spike from the rank 1 update to w. */
00060     w += v[n-1] * u;
00061 
00062     /* eliminate the spike. */
00063     *sing = false;
00064     for (j = 0; j < n-1; ++j) {
00065         if (w[j] != 0.) {
00066             /* determine a givens rotation which eliminates the */
00067             /* j-th element of the spike. */
00068             givens.makeGivens(-s(j,j), w[j]);
00069 
00070             /* apply the transformation to s and reduce the spike in w. */
00071             for (i = j; i < m; ++i) {
00072                 temp = givens.c() * s(j,i) + givens.s() * w[i];
00073                 w[i] = -givens.s() * s(j,i) + givens.c() * w[i];
00074                 s(j,i) = temp;
00075             }
00076 
00077             /* store the information necessary to recover the */
00078             /* givens rotation. */
00079             w_givens[j] = givens;
00080         } else
00081             v_givens[j] = IdentityRotation;
00082 
00083         /* test for zero diagonal elements in the output s. */
00084         if (s(j,j) == 0.) {
00085             *sing = true;
00086         }
00087     }
00088     /* move w back into the last column of the output s. */
00089     s(n-1,n-1) = w[n-1];
00090 
00091     if (s(j,j) == 0.) {
00092         *sing = true;
00093     }
00094     return;
00095 }
00096 
00097 } // end namespace internal
00098 
00099 } // end namespace Eigen


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Thu Aug 27 2015 11:59:51