Home | History | Annotate | Download | only in NonLinearOptimization 
      1 namespace Eigen {
      2 
      3 namespace internal {
      4 
      5 template <typename Scalar>
      6 void r1updt(
      7         Matrix< Scalar, Dynamic, Dynamic > &s,
      8         const Matrix< Scalar, Dynamic, 1> &u,
      9         std::vector<JacobiRotation<Scalar> > &v_givens,
     10         std::vector<JacobiRotation<Scalar> > &w_givens,
     11         Matrix< Scalar, Dynamic, 1> &v,
     12         Matrix< Scalar, Dynamic, 1> &w,
     13         bool *sing)
     14 {
     15     typedef DenseIndex Index;
     16     const JacobiRotation<Scalar> IdentityRotation = JacobiRotation<Scalar>(1,0);
     17 
     18     /* Local variables */
     19     const Index m = s.rows();
     20     const Index n = s.cols();
     21     Index i, j=1;
     22     Scalar temp;
     23     JacobiRotation<Scalar> givens;
     24 
     25     // r1updt had a broader usecase, but we dont use it here. And, more
     26     // importantly, we can not test it.
     27     assert(m==n);
     28     assert(u.size()==m);
     29     assert(v.size()==n);
     30     assert(w.size()==n);
     31 
     32     /* move the nontrivial part of the last column of s into w. */
     33     w[n-1] = s(n-1,n-1);
     34 
     35     /* rotate the vector v into a multiple of the n-th unit vector */
     36     /* in such a way that a spike is introduced into w. */
     37     for (j=n-2; j>=0; --j) {
     38         w[j] = 0.;
     39         if (v[j] != 0.) {
     40             /* determine a givens rotation which eliminates the */
     41             /* j-th element of v. */
     42             givens.makeGivens(-v[n-1], v[j]);
     43 
     44             /* apply the transformation to v and store the information */
     45             /* necessary to recover the givens rotation. */
     46             v[n-1] = givens.s() * v[j] + givens.c() * v[n-1];
     47             v_givens[j] = givens;
     48 
     49             /* apply the transformation to s and extend the spike in w. */
     50             for (i = j; i < m; ++i) {
     51                 temp = givens.c() * s(j,i) - givens.s() * w[i];
     52                 w[i] = givens.s() * s(j,i) + givens.c() * w[i];
     53                 s(j,i) = temp;
     54             }
     55         } else
     56             v_givens[j] = IdentityRotation;
     57     }
     58 
     59     /* add the spike from the rank 1 update to w. */
     60     w += v[n-1] * u;
     61 
     62     /* eliminate the spike. */
     63     *sing = false;
     64     for (j = 0; j < n-1; ++j) {
     65         if (w[j] != 0.) {
     66             /* determine a givens rotation which eliminates the */
     67             /* j-th element of the spike. */
     68             givens.makeGivens(-s(j,j), w[j]);
     69 
     70             /* apply the transformation to s and reduce the spike in w. */
     71             for (i = j; i < m; ++i) {
     72                 temp = givens.c() * s(j,i) + givens.s() * w[i];
     73                 w[i] = -givens.s() * s(j,i) + givens.c() * w[i];
     74                 s(j,i) = temp;
     75             }
     76 
     77             /* store the information necessary to recover the */
     78             /* givens rotation. */
     79             w_givens[j] = givens;
     80         } else
     81             v_givens[j] = IdentityRotation;
     82 
     83         /* test for zero diagonal elements in the output s. */
     84         if (s(j,j) == 0.) {
     85             *sing = true;
     86         }
     87     }
     88     /* move w back into the last column of the output s. */
     89     s(n-1,n-1) = w[n-1];
     90 
     91     if (s(j,j) == 0.) {
     92         *sing = true;
     93     }
     94     return;
     95 }
     96 
     97 } // end namespace internal
     98 
     99 } // end namespace Eigen
    100