eigen/unsupported/Eigen/src/NonLinearOptimization/rwupdt.h

template <typename Scalar>
void ei_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();
    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;
        }

        if (rowj == 0.)
        {
          givens[j] = JacobiRotation<Scalar>(1,0);
          continue;
        }

        /* determine a givens rotation which eliminates w(j). */
        givens[j].makeGivens(-r(j,j), rowj);

        /* 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;
    }
}